absorptioncorrutils.py 27.5 KB
Newer Older
1
2
# Mantid Repository : https://github.com/mantidproject/mantid
#
3
# Copyright © 2020 ISIS Rutherford Appleton Laboratory UKRI,
4
5
6
#   NScD Oak Ridge National Laboratory, European Spallation Source,
#   Institut Laue - Langevin & CSNS, Institute of High Energy Physics, CAS
# SPDX - License - Identifier: GPL - 3.0 +
Kendrick, Coleman's avatar
Kendrick, Coleman committed
7
from mantid.api import AnalysisDataService, WorkspaceFactory
8
from mantid.kernel import Logger, Property, PropertyManager
Zhang, Chen's avatar
Zhang, Chen committed
9
10
11
12
from mantid.simpleapi import (AbsorptionCorrection, DeleteWorkspace, Divide, Load, Multiply,
                              PaalmanPingsAbsorptionCorrection, PreprocessDetectorsToMD,
                              RenameWorkspace, SetSample, SaveNexusProcessed, UnGroupWorkspace, mtd)
import mantid.simpleapi
13
14
import numpy as np
import os
Zhang, Chen's avatar
Zhang, Chen committed
15
from functools import wraps
16

17
VAN_SAMPLE_DENSITY = 0.0721
18
19
20
_EXTENSIONS_NXS = ["_event.nxs", ".nxs.h5"]


21
22
23
# ---------------------------- #
# ----- Helper functions ----- #
# ---------------------------- #
24
def _getBasename(filename):
25
26
27
    """
    Helper function to get the filename without the path or extension
    """
28
29
30
31
32
33
34
35
    if type(filename) == list:
        filename = filename[0]
    name = os.path.split(filename)[-1]
    for extension in _EXTENSIONS_NXS:
        name = name.replace(extension, '')
    return name


Zhang, Chen's avatar
Zhang, Chen committed
36
<<<<<<< HEAD
Zhang, Chen's avatar
Zhang, Chen committed
37
<<<<<<< HEAD
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
def _getInstrName(filename, wksp=None):
    """
    Infers the instrument name from the given filename, uses wksp to fallback on
    if the instrument from the filename was invalid. In the worst case, get the
    instrument from the ConfigService.

    :param filename: Filename to use when finding instrument (ex: PG3_123 should return PG3)
    :param wksp: Optional workspace to get the instrument from
    :return: instrument name (shortname)
    """
    # First, strip off path and extensions
    base = _getBasename(filename)
    # Assume files are named as "<instr>_<id>"
    name = base.split("_")[0]
    # Check whether it is a valid instrument
    try:
        instr = mantid.kernel.ConfigService.getInstrument(name)
    except RuntimeError:
        if wksp is not None:
            # Use config service to lookup InstrumentInfo obj from workspace Instrument
            instr = mantid.kernel.ConfigService.getInstrument(wksp.getInstrument().getName())
        else:
            # Use default instrument name
            instr = mantid.kernel.ConfigService.getInstrument()
    return instr.shortName()


def _getCacheName(wkspname, wksp_prefix, cache_dir, abs_method):
66
67
68
69
    """
    Generate a MDF5 string based on given key properties.

    :param wkspname: donor workspace containing absorption correction info
70
    :param wksp_prefix: prefix to add to wkspname for caching
71
    :param cache_dir: location to store the cached absorption correction
72

73
    return fileName(full path), sha1
74
    """
Zhang, Chen's avatar
Zhang, Chen committed
75

76
    # parse algorithm used for absorption calculation
77
    # NOTE: the actual algorithm used depends on selected abs_method, therefore
78
79
80
81
82
83
84
85
    #       we are embedding the algorithm name into the SHA1 so that we can
    #       differentiate them for caching purpose
    alg_name = {
        "SampleOnly": "AbsorptionCorrection",
        "SampleAndContainer": "AbsorptionCorrection",
        "FullPaalmanPings": "PaalmanPingsAbsorptionCorrection",
    }[abs_method]

86
    # key property to generate the HASH
87
88
89
    if not mtd.doesExist(wkspname):
        raise RuntimeError(
            "Could not find workspace '{}' in ADS to create cache name".format(wkspname))
90
91
92
93
    ws = mtd[wkspname]
    # NOTE:
    #  - The query for height is tied to a beamline, which is not a good design as it
    #    will break for other beamlines
Zhang, Chen's avatar
Zhang, Chen committed
94
    property_string = [
95
96
97
98
99
100
101
102
103
104
105
106
107
        f"{key}={val}" for key, val in {
            'wavelength_min': ws.readX(0).min(),
            'wavelength_max': ws.readX(0).max(),
            "num_wl_bins": len(ws.readX(0)) - 1,
            "sample_formula": ws.run()['SampleFormula'].lastValue().strip(),
            "mass_density": ws.run()['SampleDensity'].lastValue(),
            "height_unit": ws.run()['BL11A:CS:ITEMS:HeightInContainerUnits'].lastValue(),
            "height": ws.run()['BL11A:CS:ITEMS:HeightInContainer'].lastValue(),
            "sample_container": ws.run()['SampleContainer'].lastValue().replace(" ", ""),
            "algorithm_used": alg_name,
        }.items()
    ]

Zhang, Chen's avatar
Zhang, Chen committed
108
    cache_path, signature = mantid.simpleapi.CreateCacheFilename(
109
        Prefix=wksp_prefix,
Zhang, Chen's avatar
Zhang, Chen committed
110
111
112
113
114
        OtherProperties=property_string,
        CacheDir=cache_dir,
    )

    return cache_path, signature
115
116


117
def _getCachedData(absName, abs_method, sha1, cache_file_name):
Zhang, Chen's avatar
Zhang, Chen committed
118
119
=======
def __get_cache_name(wksp_name="", cache_dir="", abs_method="SampleAndContainer"):
Zhang, Chen's avatar
Zhang, Chen committed
120
121
122
=======
def __get_cache_name(meta_wksp_name, abs_method, cache_dir=""):
>>>>>>> switch abs va call
Zhang, Chen's avatar
Zhang, Chen committed
123
124
    """generate cachefile name (full path) and sha1

Zhang, Chen's avatar
Zhang, Chen committed
125
    :param meta_wksp_name: name of workspace contains relevant meta data for hashing
Zhang, Chen's avatar
Zhang, Chen committed
126
127
128
129
130
131
    :param cache_dir: cache directory to scan/load cache data
    :param abs_method: method used to perform the absorption calculation
    
    return cachefile_name: full path of the cache file
           sha1: MD5 value based on selected property
>>>>>>> try the decorator approach
132
    """
Zhang, Chen's avatar
Zhang, Chen committed
133
    # grab the workspace
Zhang, Chen's avatar
Zhang, Chen committed
134
135
    if meta_wksp_name in mtd:
        ws = mtd[meta_wksp_name]
Zhang, Chen's avatar
Zhang, Chen committed
136
137
    else:
        raise ValueError(
Zhang, Chen's avatar
Zhang, Chen committed
138
            f"Cannot find workspace {meta_wksp_name} to extract meta data for hashing, aborting")
139

Zhang, Chen's avatar
Zhang, Chen committed
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
    # requires cache_dir
    if cache_dir == "":
        cache_filename, signature = "", ""
    else:
        # generate the property string for hashing
        property_string = [
            f"{key}={val}" for key, val in {
                'wavelength_min': ws.readX(0).min(),
                'wavelength_max': ws.readX(0).max(),
                "num_wl_bins": len(ws.readX(0)) - 1,
                "sample_formula": ws.run()['SampleFormula'].lastValue().strip(),
                "mass_density": ws.run()['SampleDensity'].lastValue(),
                "height_unit": ws.run()['BL11A:CS:ITEMS:HeightInContainerUnits'].lastValue(),
                "height": ws.run()['BL11A:CS:ITEMS:HeightInContainer'].lastValue(),
                "sample_container": ws.run()['SampleContainer'].lastValue().replace(" ", ""),
                "abs_method": abs_method,
            }.items()
        ]

        # use mantid build-in alg to generate the cache filename and sha1
        cache_filename, signature = mantid.simpleapi.CreateCacheFilename(
            OtherProperties=property_string,
            CacheDir=cache_dir,
        )

    return cache_filename, signature


def __load_cached_data(cache_file_name, sha1, abs_method=""):
    """try to load cached data from memory and disk

    :param abs_method: absorption calculation method
    :param sha1: SHA1 that identify cached workspace
    :param cache_file_name: cache file name to search

    return abs_wksp_sample,
           abs_wksp_container
177
    """
Zhang, Chen's avatar
Zhang, Chen committed
178
179
180
    # init
    abs_wksp_sample, abs_wksp_container = "", ""
    found_abs_wksp_sample, found_abs_wksp_container = False, False
181

Zhang, Chen's avatar
Zhang, Chen committed
182
183
184
    # step_0: depending on the abs_method, prefix will be different
    # NOTE:
    #  This is our internal naming convention for abs workspaces
185
    if abs_method == "SampleOnly":
Zhang, Chen's avatar
Zhang, Chen committed
186
187
        abs_wksp_sample = f"abs_ass_{sha1}"
        found_abs_wksp_container = True
188
    elif abs_method == "SampleAndContainer":
Zhang, Chen's avatar
Zhang, Chen committed
189
190
        abs_wksp_sample = f"abs_ass_{sha1}"
        abs_wksp_container = f"abs_acc_{sha1}"
191
    elif abs_method == "FullPaalmanPings":
Zhang, Chen's avatar
Zhang, Chen committed
192
193
        abs_wksp_sample = f"abs_assc_{sha1}"
        abs_wksp_container = f"abs_ac_{sha1}"
194
    else:
195
        raise ValueError("Unrecognized absorption correction method '{}'".format(abs_method))
196

Zhang, Chen's avatar
Zhang, Chen committed
197
198
199
200
201
202
    # step_1: check memory
    found_abs_wksp_sample = mtd.doesExist(abs_wksp_sample)
    found_abs_wksp_container = mtd.doesExist(abs_wksp_container)

    # step_2: load from disk if either is not found in memory
    if (not found_abs_wksp_sample) or (not found_abs_wksp_container):
Zhang, Chen's avatar
Zhang, Chen committed
203
        if os.path.exists(cache_file_name):
Zhang, Chen's avatar
Zhang, Chen committed
204
            wsntmp = f"tmpwsg"
Zhang, Chen's avatar
Zhang, Chen committed
205
206
            Load(Filename=cache_file_name, OutputWorkspace=wsntmp)
            wstype = mtd[wsntmp].id()
207
            if wstype == "Workspace2D":
Zhang, Chen's avatar
Zhang, Chen committed
208
                RenameWorkspace(InputWorkspace=wsntmp, OutputWorkspace=abs_wksp_sample)
Zhang, Chen's avatar
Zhang, Chen committed
209
210
211
            elif wstype == "WorkspaceGroup":
                UnGroupWorkspace(InputWorkspace=wsntmp)
            else:
212
                raise ValueError(f"Unsupported cached workspace type: {wstype}")
213

Zhang, Chen's avatar
Zhang, Chen committed
214
215
216
    # step_3: check memory again
    found_abs_wksp_sample = mtd.doesExist(abs_wksp_sample)
    found_abs_wksp_container = mtd.doesExist(abs_wksp_container)
217

Zhang, Chen's avatar
Zhang, Chen committed
218
219
    abs_wksp_sample = abs_wksp_sample if found_abs_wksp_sample else ""
    abs_wksp_container = abs_wksp_container if found_abs_wksp_container else ""
220

Zhang, Chen's avatar
Zhang, Chen committed
221
222
223
224
225
226
227
    return abs_wksp_sample, abs_wksp_container


# NOTE:
#  In order to use the decorator, we must have consistent naming
#  or kwargs as this is probably the most reliable way to get
#  the desired data piped in multiple location
Zhang, Chen's avatar
Zhang, Chen committed
228
229
#  -- bare minimum signaure of the function
#    func(wksp_name: str, abs_method:str, cache_dir="")
Zhang, Chen's avatar
Zhang, Chen committed
230
231
232
233
def abs_cache(func):
    """decorator to make the caching process easier"""
    @wraps(func)
    def inner(*args, **kwargs):
Zhang, Chen's avatar
Zhang, Chen committed
234
235
236
237
238
        # unpack key arguments
        wksp_name = args[0]
        abs_method = args[1]
        cache_dir = kwargs.get("cache_dir", "")

Zhang, Chen's avatar
Zhang, Chen committed
239
        # prompt return if no cache_dir specified
Zhang, Chen's avatar
Zhang, Chen committed
240
        if cache_dir == "":
Zhang, Chen's avatar
Zhang, Chen committed
241
242
243
244
            return func(*args, **kwargs)

        # step_1: generate the SHA1 and cachefile name
        #         baseon given kwargs
Zhang, Chen's avatar
Zhang, Chen committed
245
        cache_filename, signature = __get_cache_name(wksp_name, abs_method, cache_dir)
Zhang, Chen's avatar
Zhang, Chen committed
246
247
248

        # step_2: try load the cached data
        abs_wksp_sample, abs_wksp_container = __load_cached_data(cache_filename, signature,
Zhang, Chen's avatar
Zhang, Chen committed
249
                                                                 abs_method)
Zhang, Chen's avatar
Zhang, Chen committed
250
251

        # step_3: calculation
Zhang, Chen's avatar
Zhang, Chen committed
252
        if (abs_method == "SampleOnly") and (abs_wksp_sample != ""):
Zhang, Chen's avatar
Zhang, Chen committed
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
            return abs_wksp_sample
        else:
            if (abs_wksp_sample != "") and (abs_wksp_container != ""):
                # cache is available in memory now, skip calculation
                return abs_wksp_sample, abs_wksp_container
            else:
                # no cache found, need calculation
                abs_wksp_sample_nosha, abs_wksp_container_nosha = func(*args, **kwargs)
                # need to add the sha in the workspace name
                mantid.simpleapi.RenameWorkspace(abs_wksp_sample_nosha, abs_wksp_sample)
                SaveNexusProcessed(InputWorkspace=abs_wksp_sample, Filename=cache_filename)

                if abs_wksp_container_nosha != "":
                    mantid.simpleapi.RenameWorkspace(abs_wksp_container_nosha, abs_wksp_container)
                    SaveNexusProcessed(InputWorkspace=abs_wksp_container,
                                       Filename=cache_filename,
                                       Append=True)
                else:
                    abs_wksp_container = ""

                return abs_wksp_sample, abs_wksp_container

    return inner
276
277
278
279
280


# ----------------------------- #
# ---- Core functionality ----- #
# ------------------------------#
281
def calculate_absorption_correction(
282
283
284
285
286
287
288
289
290
291
292
293
        filename,
        abs_method,
        props,
        sample_formula,
        mass_density,
        number_density=Property.EMPTY_DBL,
        container_shape="PAC06",
        num_wl_bins=1000,
        element_size=1,
        metaws=None,
        cache_dir="",
        prefix="FILENAME",
294
):
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
    """The absorption correction is applied by (I_s - I_c*k*A_csc/A_cc)/A_ssc for pull Paalman-Ping

    If no cross-term then I_s/A_ss - I_c/A_cc

    Therefore this will return 2 workspace, one for correcting the
    sample (A_s) and one for the container (A_c) depending on the
    absorption method, that will be passed to _focusAndSum and
    therefore AlignAndFocusPowderFromFiles.

    If SampleOnly then

    A_s = A_ss
    A_c = None

    If SampleAndContainer then

    A_s = A_ss
    A_c = A_cc

    If FullPaalmanPings then
    A_s = A_ssc
    A_c = A_cc*A_ssc/A_csc

    This will then return (A_s, A_c)
319
320
321

    :param filename: File to be used for absorption correction
    :param abs_method: Type of absorption correction: None, SampleOnly, SampleAndContainer, FullPaalmanPings
322
    :param props: PropertyManager of run characterizations, obtained from PDDetermineCharacterizations
323
324
325
326
327
328
    :param sample_formula: Sample formula to specify the Material for absorption correction
    :param mass_density: Mass density of the sample to specify the Material for absorption correction
    :param number_density: Optional number density of sample to be added to the Material for absorption correction
    :param container_shape: Shape definition of container, such as PAC06.
    :param num_wl_bins: Number of bins for calculating wavelength
    :param element_size: Size of one side of the integration element cube in mm
329
    :param metaws: Optional workspace containing metadata to use instead of reading from filename
330
    :param cache_dir: cache directory for storing cached absorption correction workspace
331
    :param prefix: How the prefix of cache file is determined - FILENAME to use file, or SHA prefix
332

333
    :return:
334
        Two workspaces (A_s, A_c) names
335
336
337
338
    """
    if abs_method == "None":
        return None, None

339
    material = {"ChemicalFormula": sample_formula, "SampleMassDensity": mass_density}
340
341
342
343
344
345

    if number_density != Property.EMPTY_DBL:
        material["SampleNumberDensity"] = number_density

    environment = {'Name': 'InAir', 'Container': container_shape}

346
347
348
349
350
351
    donorWS = create_absorption_input(filename,
                                      props,
                                      num_wl_bins,
                                      material=material,
                                      environment=environment,
                                      metaws=metaws)
352

Zhang, Chen's avatar
Zhang, Chen committed
353
    absName = '{}_abs_correction'.format(_getBasename(filename))
354

355
    if cache_dir == "":
356
357
358
359
        # no caching activity if no cache directory is provided
        return calc_absorption_corr_using_wksp(donorWS, abs_method, element_size, absName)
    else:
        # -- Caching -- #
360
361
362
363
364
        # -- Generate cache file prefix: defaults to use base filename unless prefix arg is set to SHA
        cache_prefix = _getInstrName(filename, donorWS)
        if prefix != "SHA":
            # Append the filename to the instrument name for the cache file
            cache_prefix = cache_prefix + _getBasename(filename).replace(cache_prefix, '')
365
        # -- Generate the cache file name based on
366
        #    - input filename (embedded in donorWS) or SHA depending on prefix determined above
367
        #    - SNSPowderReduction Options (mostly passed in as args)
368
        cache_filename, sha1 = _getCacheName(donorWS, cache_prefix, cache_dir, abs_method)
369
        log.information(f"SHA1: {sha1}")
370
371
372
373
374
375
        # -- Try to use cache
        #    - if cache is found, wsn_as and wsn_ac will be valid string (workspace name)
        #      - already exist in memory
        #      - load from cache nxs file
        #    - if cache is not found, wsn_as and wsn_ac will both be None
        #      - standard calculation will be kicked off as before
376
377
378
379
380
381
        # Update absName with the cache prefix to find workspaces in memory
        if prefix == "SHA":
            absName = cache_prefix + "_" + sha1 + "_abs_correction"
        else:
            absName = cache_prefix + "_abs_correction"
        wsn_as, wsn_ac = _getCachedData(absName, abs_method, sha1, cache_filename)
382
383

        # NOTE:
384
385
386
387
        # -- one algorithm with three very different behavior, why not split them to
        #    to make the code cleaner, also the current design will most likely leads
        #    to severe headache down the line
        log.information(f"For current analysis using {abs_method}")
388
        if (abs_method == "SampleOnly") and (wsn_as != ""):
389
390
            # first deal with special case where we only care about the sample absorption
            log.information(f"-- Located cached workspace, {wsn_as}")
391
            # NOTE:
392
393
394
            #  Nothing to do here, since
            #  - wsn_as is already loaded by _getCachedData
            #  - wsn_ac is already set to None by _getCachedData.
395
        else:
396
            if (wsn_as == "") or (wsn_ac == ""):
397
398
                log.information(f"-- Cannot locate all necessary cache, start from scrach")
                wsn_as, wsn_ac = calc_absorption_corr_using_wksp(donorWS, abs_method, element_size,
399
                                                                 absName)
400
401
402
403
404
405
406
                # NOTE:
                #  We need to set the SHA1 first, then save.
                #  Because the final check is always comparing SHA1 of given
                #  workspace.
                # set the SHA1 to workspace in memory (for in-memory cache search)
                mtd[wsn_as].mutableRun()["absSHA1"] = sha1
                # case SampleOnly is the annoying one
407
                if wsn_ac is not None:
408
                    mtd[wsn_ac].mutableRun()["absSHA1"] = sha1
409

410
411
412
                # save the cache to file (for hard-disk cache)
                SaveNexusProcessed(InputWorkspace=wsn_as, Filename=cache_filename)
                # case SampleOnly is the annoying one
413
                if wsn_ac is not None:
414
415
416
417
                    SaveNexusProcessed(InputWorkspace=wsn_ac, Filename=cache_filename, Append=True)
            else:
                # found the cache, let's use the cache instead
                log.information(f"-- Locate cached sample absorption correction: {wsn_as}")
418
                log.information(f"-- Locate cached container absorption correction: {wsn_ac}")
419
420

    return wsn_as, wsn_ac
Zhang, Chen's avatar
Zhang, Chen committed
421
422
423
424
425
426
427
428
429
430
431
    return calc_absorption_corr_using_wksp(donorWS,
                                           abs_method,
                                           element_size,
                                           absName,
                                           cache_dir=cache_dir)


@abs_cache
def calc_absorption_corr_using_wksp(donor_wksp,
                                    abs_method,
                                    element_size=1,
Zhang, Chen's avatar
Zhang, Chen committed
432
                                    prefix_name="",
Zhang, Chen's avatar
Zhang, Chen committed
433
                                    cache_dir=""):
434
435
436
437
438
439
440
441
    """
    Calculates absorption correction on the specified donor workspace. See the documentation
    for the ``calculate_absorption_correction`` function above for more details.

    :param donor_wksp: Input workspace to compute absorption correction on
    :param abs_method: Type of absorption correction: None, SampleOnly, SampleAndContainer, FullPaalmanPings
    :param element_size: Size of one side of the integration element cube in mm
    :param prefix_name: Optional prefix of the output workspaces, default is the donor_wksp name.
Zhang, Chen's avatar
Zhang, Chen committed
442
443
    :param cache_dir: Cache directory to store cached abs workspace.

444
445
    :return: Two workspaces (A_s, A_c), the first for the sample and the second for the container
    """
Zhang, Chen's avatar
Zhang, Chen committed
446
447
448
    log = Logger('calc_absorption_corr_using_wksp')
    if cache_dir != "":
        log.information(f"Storing cached data in {cache_dir}")
449
450

    if abs_method == "None":
Zhang, Chen's avatar
Zhang, Chen committed
451
        return "", ""
452
453
454
455
456
457
458
459
460

    if isinstance(donor_wksp, str):
        if not mtd.doesExist(donor_wksp):
            raise RuntimeError("Specified donor workspace not found in the ADS")
        donor_wksp = mtd[donor_wksp]

    absName = donor_wksp.name()
    if prefix_name != '':
        absName = prefix_name
461
462

    if abs_method == "SampleOnly":
463
        AbsorptionCorrection(donor_wksp,
464
465
466
                             OutputWorkspace=absName + '_ass',
                             ScatterFrom='Sample',
                             ElementSize=element_size)
Zhang, Chen's avatar
Zhang, Chen committed
467
        return absName + '_ass', ""
468
    elif abs_method == "SampleAndContainer":
469
        AbsorptionCorrection(donor_wksp,
470
471
472
                             OutputWorkspace=absName + '_ass',
                             ScatterFrom='Sample',
                             ElementSize=element_size)
473
        AbsorptionCorrection(donor_wksp,
474
475
476
477
                             OutputWorkspace=absName + '_acc',
                             ScatterFrom='Container',
                             ElementSize=element_size)
        return absName + '_ass', absName + '_acc'
478
    elif abs_method == "FullPaalmanPings":
479
        PaalmanPingsAbsorptionCorrection(donor_wksp,
480
481
482
483
484
485
486
487
488
                                         OutputWorkspace=absName,
                                         ElementSize=element_size)
        Multiply(LHSWorkspace=absName + '_acc',
                 RHSWorkspace=absName + '_assc',
                 OutputWorkspace=absName + '_ac')
        Divide(LHSWorkspace=absName + '_ac',
               RHSWorkspace=absName + '_acsc',
               OutputWorkspace=absName + '_ac')
        return absName + '_assc', absName + '_ac'
489
    else:
490
        raise ValueError("Unrecognized absorption correction method '{}'".format(abs_method))
491
492


493
def create_absorption_input(
494
495
496
497
498
499
500
501
502
        filename,
        props,
        num_wl_bins=1000,
        material=None,
        geometry=None,
        environment=None,
        opt_wl_min=0,
        opt_wl_max=Property.EMPTY_DBL,
        metaws=None,
503
):
504
505
    """
    Create an input workspace for carpenter or other absorption corrections
506
507

    :param filename: Input file to retrieve properties from the sample log
508
    :param props: PropertyManager of run characterizations, obtained from PDDetermineCharacterizations
509
510
511
512
    :param num_wl_bins: The number of wavelength bins used for absorption correction
    :param material: Optional material to use in SetSample
    :param geometry: Optional geometry to use in SetSample
    :param environment: Optional environment to use in SetSample
513
514
    :param opt_wl_min: Optional minimum wavelength. If specified, this is used instead of from the props
    :param opt_wl_max: Optional maximum wavelength. If specified, this is used instead of from the props
515
    :param metaws: Optional workspace name with metadata to use for donor workspace instead of reading from filename
516
    :return: Name of the donor workspace created
517
    """
518
519
520
521
522
    if props is None:
        raise RuntimeError("props is required to create donor workspace, props is None")

    if not isinstance(props, PropertyManager):
        raise RuntimeError("props must be a PropertyManager object")
523
524
525

    log = Logger('CreateAbsorptionInput')

526
527
528
529
    # Load from file if no workspace with metadata has been given, otherwise avoid a duplicate load with the metaws
    absName = metaws
    if metaws is None:
        absName = '__{}_abs'.format(_getBasename(filename))
Zhang, Chen's avatar
Zhang, Chen committed
530
531
532
533
534
        allowed_log = " ".join([
            'SampleFormula', 'SampleDensity', "BL11A:CS:ITEMS:HeightInContainerUnits",
            "SampleContainer"
        ])
        Load(Filename=filename, OutputWorkspace=absName, MetaDataOnly=True, AllowList=allowed_log)
535
536

    # first attempt to get the wavelength range from the properties file
537
    wl_min, wl_max = props['wavelength_min'].value, props['wavelength_max'].value
538
539
540
541
542
543
544
545
    # override that with what was given as parameters to the algorithm
    if opt_wl_min > 0.:
        wl_min = opt_wl_min
    if opt_wl_max != Property.EMPTY_DBL:
        wl_max = opt_wl_max

    # if it isn't found by this point, guess it from the time-of-flight range
    if (wl_min == wl_max == 0.):
546
547
        tof_min = props['tof_min'].value
        tof_max = props['tof_max'].value
548
549
550
551
552
553
554
        if tof_min >= 0. and tof_max > tof_min:
            log.information('TOF range is {} to {} microseconds'.format(tof_min, tof_max))

            # determine L1
            instr = mtd[absName].getInstrument()
            L1 = instr.getSource().getDistance(instr.getSample())
            # determine L2 range
555
556
557
            PreprocessDetectorsToMD(InputWorkspace=absName,
                                    OutputWorkspace=absName + '_dets',
                                    GetMaskState=False)
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
            L2 = mtd[absName + '_dets'].column('L2')
            Lmin = np.min(L2) + L1
            Lmax = np.max(L2) + L1
            DeleteWorkspace(Workspace=absName + '_dets')

            log.information('Distance range is {} to {} meters'.format(Lmin, Lmax))

            # wavelength is h*TOF / m_n * L  values copied from Kernel/PhysicalConstants.h
            usec_to_sec = 1.e-6
            meter_to_angstrom = 1.e10
            h_m_n = meter_to_angstrom * usec_to_sec * 6.62606896e-34 / 1.674927211e-27
            wl_min = h_m_n * tof_min / Lmax
            wl_max = h_m_n * tof_max / Lmin

    # there isn't a good way to guess it so error out
    if wl_max <= wl_min:
        DeleteWorkspace(Workspace=absName)  # no longer needed
        raise RuntimeError('Invalid wavelength range min={}A max={}A'.format(wl_min, wl_max))
    log.information('Using wavelength range min={}A max={}A'.format(wl_min, wl_max))

578
579
580
581
    absorptionWS = WorkspaceFactory.create(mtd[absName],
                                           NVectors=mtd[absName].getNumberHistograms(),
                                           XLength=num_wl_bins + 1,
                                           YLength=num_wl_bins)
582
583
584
585
586
587
588
589
590
591
592
593
    xaxis = np.arange(0., float(num_wl_bins + 1)) * (wl_max - wl_min) / (num_wl_bins) + wl_min
    for i in range(absorptionWS.getNumberHistograms()):
        absorptionWS.setX(i, xaxis)
    absorptionWS.getAxis(0).setUnit('Wavelength')

    # this effectively deletes the metadata only workspace
    AnalysisDataService.addOrReplace(absName, absorptionWS)

    # Set ChemicalFormula, and either SampleMassDensity or Mass, if SampleMassDensity not set
    if material is not None:
        if (not material['ChemicalFormula']) and ("SampleFormula" in absorptionWS.run()):
            material['ChemicalFormula'] = absorptionWS.run()['SampleFormula'].lastValue().strip()
594
595
596
597
        if ("SampleMassDensity" not in material
                or not material['SampleMassDensity']) and ("SampleDensity" in absorptionWS.run()):
            if (absorptionWS.run()['SampleDensity'].lastValue() !=
                    1.0) and (absorptionWS.run()['SampleDensity'].lastValue() != 0.0):
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
                material['SampleMassDensity'] = absorptionWS.run()['SampleDensity'].lastValue()
            else:
                material['Mass'] = absorptionWS.run()['SampleMass'].lastValue()

    # Set height for computing density if height not set
    if geometry is None:
        geometry = {}

    if geometry is not None:
        if "Height" not in geometry or not geometry['Height']:
            # Check units - SetSample expects cm
            if absorptionWS.run()['BL11A:CS:ITEMS:HeightInContainerUnits'].lastValue() == "mm":
                conversion = 0.1
            elif absorptionWS.run()['BL11A:CS:ITEMS:HeightInContainerUnits'].lastValue() == "cm":
                conversion = 1.0
            else:
614
615
616
                raise ValueError(
                    "HeightInContainerUnits expects cm or mm; specified units not recognized: ",
                    absorptionWS.run()['BL11A:CS:ITEMS:HeightInContainerUnits'].lastValue())
617

618
619
            geometry['Height'] = absorptionWS.run()['BL11A:CS:ITEMS:HeightInContainer'].lastValue(
            ) * conversion
620
621
622
623

    # Set container if not set
    if environment is not None:
        if environment['Container'] == "":
624
625
            environment['Container'] = absorptionWS.run()['SampleContainer'].lastValue().replace(
                " ", "")
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648

    # Make sure one is set before calling SetSample
    if material or geometry or environment is not None:
        setup_sample(absName, material, geometry, environment)

    return absName


def setup_sample(donor_ws, material, geometry, environment):
    """
    Calls SetSample with the associated sample and container material and geometry for use
    in creating an input workspace for an Absorption Correction algorithm
    :param donor_ws:
    :param material:
    :param geometry:
    :param environment:
    """

    # Set the material, geometry, and container info
    SetSample(InputWorkspace=donor_ws,
              Material=material,
              Geometry=geometry,
              Environment=environment)