Update POLDI (PSI) auto-correlation algorithm for new 2D detector

[RichardWaiteSTFC]

Is your feature request related to a problem? Please describe.
POLDI have a new 2D detector, but the old auto-correlation algorithm POLDIAutoCorellation does not work for this type of detector (assumes all pixels at two-theta = 90 deg). Also it is not using the pixel positions in the IDF for the unit conversions.

Describe the solution you'd like
New python (TBC) algorithm that supports auto-correlation method that uses detector positions as in the IDF

Prototype code I wrote in python

# import mantid algorithms, numpy and matplotlib
from mantid.simpleapi import *
import matplotlib.pyplot as plt
import numpy as np
from mantid.kernel import UnitConversion, DeltaEModeType, UnitParams, UnitParametersMap
from os import path
from time import time

# NEW POLDI
###########

data_dir = r"C:\Users\xhg73778\Downloads\ForRichardPOLDI"
# data_dir = "/babylon/Public/RWaite/POLDI/"

ws = LoadEmptyInstrument(Filename=path.join(data_dir, "POLDI_Definition_newDet_10binsPchannel - kopie.xml")) #  - kopie
d = np.loadtxt(path.join(data_dir, "summedDiffractograms4480x500", "summedDiffractograms.txt"))
xmax = 3000 # 1497
bin_width = xmax/d.shape[-1]
ws = Rebin(ws, Params=f"0,{bin_width},{xmax}", OutputWorkspace=ws.name())
for ispec in range(ws.getNumberHistograms()):
    ws.setY(ispec, d[ispec,:])
    
# add some logs (will eventually be set in file)
AddSampleLog(Workspace=ws, LogName="chopperspeed", LogText="5000.0", LogType="Number")
AddSampleLog(Workspace=ws, LogName="ChopperPhase", LogText="89.8", LogType="Number")
# set t0 (would normally live in parameter file)
SetInstrumentParameter(Workspace=ws, ComponentName='chopper', ParameterName='t0', 
                       ParameterType='Number', Value="0.025701")
SetInstrumentParameter(Workspace=ws, ComponentName='chopper', ParameterName='t0_const', 
                       ParameterType='Number', Value="85.0")

# get instrument/chopper settings
#################################

inst = ws.getInstrument()
source = inst.getSource()  # might not need these
sample = inst.getSample()
chopper = inst.getComponentByName("chopper")
l_chop = (chopper.getPos() - source.getPos()).norm()
t0 = chopper.getNumberParameter("t0")[0] 
t0_const = chopper.getNumberParameter("t0_const")[0]
chopper_speed = ws.run().getPropertyAsSingleValue("chopperspeed")  # rpm
chopper_phase = ws.run().getPropertyAsSingleValue("ChopperPhase")  # Doesn't appear to be used?
# chopper has 32 slits, 8 slits per section (4*8 - quarter repeated 4 times).
cycle_time = 60.0 / (4.0 * chopper_speed) * 1.0e6; # mus ?
zero_offset = t0*cycle_time + t0_const # opening time of first slit?
# get chopper offsetts in time from start of a given pulse
nslits = chopper.nelements()
slit_offsets = np.array([chopper[islit].getPos()[0]*cycle_time + t0 for islit in range(nslits)]) # offset in time

# get d-spacing array
lambda_min, lambda_max = 1.1, 5
si = ws.spectrumInfo()
tths = np.array([si.twoTheta(ispec) for ispec in range(ws.getNumberHistograms())])
tth_min, tth_max = tths.min(), tths.max()
dspac_min = lambda_min/(2*np.sin(tth_max/2))
dspac_max = lambda_max/(2*np.sin(tth_min/2))
nbins_dspac = 2850 # how to calculate this? Think this is meant to be 5531
dspacs = np.linspace(dspac_min, dspac_max, nbins_dspac) 

# get npulses to include in calc
start = time()
l2s = [si.l2(ispec) for ispec in range(ws.getNumberHistograms())]
imax = np.argmax(l2s)
params = UnitParametersMap()
params[UnitParams.l2] = l2s[imax]
params[UnitParams.twoTheta] = tths[imax]
# calc time of longest wavelength neutron to reach detector with largest total path length from first slit opening
time_max = UnitConversion.run("Wavelength", "TOF", lambda_max, si.l1() - l_chop, DeltaEModeType.Elastic, params) + slit_offsets[0]
npulses = int(time_max // cycle_time)  # 14

# evaluate intermediate corellation spectrum (Eq. 7 in POLDI concept paper)
ipulses = np.arange(npulses)[:,None]
offsets = (ipulses*cycle_time + slit_offsets).flatten()
inter_corr = np.zeros((len(dspacs), len(offsets)), dtype=float) # npulses*nslits
# sum over all spectra
for ispec in range(ws.getNumberHistograms()):
    print("ispec = ", ispec, ispec/ws.getNumberHistograms())
    params[UnitParams.l2] = l2s[ispec]
    params[UnitParams.twoTheta] = si.twoTheta(ispec)
    tof_d1Ang = UnitConversion.run("dSpacing", "TOF", 1.0, si.l1() - l_chop, DeltaEModeType.Elastic, params)
    itimes = np.mod(tof_d1Ang*dspacs[:,None] + offsets + t0_const, cycle_time) / bin_width
    ibins = np.floor(itimes).astype(int)
    ibins_plus = ibins + 1
    ibins_plus[ibins_plus > ws.blocksize() - 1 ] = 0
    y = ws.readY(ispec)
    inter_corr += (ibins + 1 - itimes)*y[ibins] + (itimes - ibins)*y[ibins_plus]

corr = 1/np.sum(1/inter_corr, axis=1)
# corr = np.mean(inter_corr, axis=1)

print("Elapsed time = ", time() - start)

fig,ax = plt.subplots()
ax.plot(dspacs, corr, '-ok')
ax.set_ylabel("Correlation Intensity (a.u.)")
ax.set_xlabel('d-spacing (Ang)')
fig.show()

fig,ax = plt.subplots()
ax.plot(2*np.pi/dspacs, corr, '-ok', markersize=3)
ax.set_ylabel("Correlation Intensity (a.u.)")
ax.set_xlabel('Q (1/Ang)')
fig.show()

Produces