SCDCalibratePanels v2¶

# necessary import
import numpy as np
from collections import namedtuple
from mantid.simpleapi import *
from mantid.geometry import CrystalStructure

# generate synthetic testing data
def convert(dictionary):
    return namedtuple('GenericDict', dictionary.keys())(**dictionary)

# lattice constant for Si
# data from Mantid web documentation
lc_silicon = {
    "a": 5.431,  # A
    "b": 5.431,  # A
    "c": 5.431,  # A
    "alpha": 90,  # deg
    "beta": 90,  # deg
    "gamma": 90,  # deg
}
silicon = convert(lc_silicon)
cs_silicon = CrystalStructure(
    f"{silicon.a} {silicon.b} {silicon.c}",
    "F d -3 m",
    "Si 0 0 0 1.0 0.05",
)

# Generate simulated workspace for CORELLI
CreateSimulationWorkspace(
    Instrument='CORELLI',
    BinParams='1,100,10000',
    UnitX='TOF',
    OutputWorkspace='cws',
)
cws = mtd['cws']

# Set the UB matrix for the sample
# u, v is the critical part, we can start with the
# ideal position
SetUB(
    Workspace="cws",
    u='1,0,0',  # vector along k_i, when goniometer is at 0
    v='0,1,0',  # in-plane vector normal to k_i, when goniometer is at 0
    **lc_silicon,
)

# set the crystal structure for virtual workspace
cws.sample().setCrystalStructure(cs_silicon)

# tweak L1
dL1 = 1.414e-2  # 1.414cm
MoveInstrumentComponent(
    Workspace='cws',
    ComponentName='moderator',
    'X'=0, 'Y'=0, 'Z'=dL1,
    RelativePosition=true,
)

# Generate predicted peak workspace
dspacings = convert({'min': 1.0, 'max': 10.0})
wavelengths = convert({'min': 0.8, 'max': 2.9})

# Collect peaks over a range of omegas
CreatePeaksWorkspace(OutputWorkspace='pws')
omegas = range(0, 180, 3)

for omega in tqdm(omegas):
    SetGoniometer(
        Workspace="cws",
        Axis0=f"{omega},0,1,0,1",
    )

    PredictPeaks(
        InputWorkspace='cws',
        WavelengthMin=wavelengths.min,
        wavelengthMax=wavelengths.max,
        MinDSpacing=dspacings.min,
        MaxDSpacing=dspacings.max,
        ReflectionCondition='All-face centred',
        OutputWorkspace='_pws',
    )

    CombinePeaksWorkspaces(
        LHSWorkspace="_pws",
        RHSWorkspace="pws",
        OutputWorkspace="pws",
    )

pws = mtd['pws']

# move the source back to make PWS forget the answer
MoveInstrumentComponent(
    Workspace='pws',
    ComponentName='moderator',
    'X'=0, 'Y'=0, 'Z'=-dL1,
    RelativePosition=true,
)

# run the calibration on pws
# similar to actual calibration, where
#   1. the peaks in side pws knows the correct L1, but info is embeded in Qsamples
#   2. the recored L1 in instrument Info is the default engineering value
SCDCalibratePanels(
    PeakWorkspace="pws",
    a=silicon.a,
    b=silicon.b,
    c=silicon.c,
    alpha=silicon.alpha,
    beta=silicon.beta,
    gamma=silicon.gamma,
    CalibrateT0=False,
    CalibrateL1=True,
    CalibrateBanks=False,
    OutputWorkspace="testCaliTable",
    XmlFilename="test.xml",
)