Consider phase offset for disk on Merlin in PyChop

docs/source/release/v6.10.0/Direct_Geometry/General/New_features/37188.rst

@@ -0,0 +1 @@
+- In PyChop, the phase offset for the disk chopper controller on MERLIN can now be modified through the yaml file for MERLIN. Furthermore, the default value has been changed from `1500` to `17000`.

qt/python/mantidqtinterfaces/mantidqtinterfaces/PyChop/PyChopGui.py

@@ -241,7 +241,7 @@ def _merlin_chopper(self):
                 self.widgets["MultiRepCheck"].setEnabled(True)
                 self.tabs.setTabEnabled(self.tdtabID, True)
                 if self.instSciAct.isChecked():
-                    self.widgets["Chopper0Phase"]["Edit"].setText("1500")
+                    self.widgets["Chopper0Phase"]["Edit"].setText("17000")
                     self.widgets["Chopper0Phase"]["Label"].setText("Disk chopper phase delay time")
                     self.widgets["Chopper0Phase"]["Edit"].show()
                     self.widgets["Chopper0Phase"]["Label"].show()
@@ -762,11 +762,11 @@ def update_script(self):
             ei = self.engine.getEi()
             out = self.engine.getWidths()
             res = out["Energy"]
-            percent = res / ei * 100
+            percent = res[0] / ei * 100
             chop_width = out["chopper"]
             mod_width = out["moderator"]
-            new_str = "\nEi is %6.2f meV, resolution is %6.2f ueV, percentage resolution is %6.3f\n" % (ei, res * 1000, percent)
-            new_str += "FWHM at sample from chopper and moderator are %6.2f us, %6.2f us\n" % (chop_width, mod_width)
+            new_str = "\nEi is %6.2f meV, resolution is %6.2f ueV, percentage resolution is %6.3f\n" % (ei, res[0] * 1000, percent)
+            new_str += "FWHM at sample from chopper and moderator are %6.2f us, %6.2f us\n" % (chop_width[0], mod_width[0])
         self.scredt.append(new_str)
 
     def onHelp(self):

qt/python/mantidqtinterfaces/test/PyChopInstrumentTest.py

@@ -14,6 +14,7 @@
 import numpy as np
 
 from pychop.Instruments import Instrument
+from pychop import MulpyRep
 
 
 class PyChopInstrumentTests(unittest.TestCase):
@@ -412,7 +413,7 @@ def test_pychop_numerics(self):
         for inst, ch, frq, ei, res0, flux0 in zip(instruments, choppers, freqs, eis, ref_res, ref_flux):
             res, flux = Instrument.calculate(inst, ch, frq, ei, 0)
             np.testing.assert_allclose(res[0], res0, rtol=1e-4, atol=0)
-            np.testing.assert_allclose(flux[0], flux0, rtol=1e-4, atol=0)
+            np.testing.assert_allclose(flux, flux0, rtol=1e-4, atol=0)
 
     def test_erange(self):
         # Tests MARI raises if Ei outside range [0, 180meV] with G chopper only ('S' chopper ok up to 1000meV)
@@ -423,6 +424,20 @@ def test_erange(self):
             obj.setEi(500)
         obj.setEi(180)
 
+    def test_phase_offset(self):
+        # Tests calcChopTimes applies new phase offset parameter
+        instpars = [[9.3, 9.995], [1, 2], None, [950, 10.0], [64, 10.0], [250, 290.0], [1, 1], 2.5, 1.925, [50, 1], 0, 0.9, [0]]
+        phaseoffset = None
+        Eis_none, _, _, _, _ = MulpyRep.calcChopTimes(50.0, [50.0, 400.0], instpars, [17000], phaseoffset)
+        self.assertAlmostEqual(Eis_none[5], 10.71281984, places=7)
+        self.assertAlmostEqual(Eis_none[6], 7.70630740, places=7)
+        self.assertAlmostEqual(Eis_none[7], 5.80834507, places=7)
+        phaseoffset = 4500
+        Eis_offset, _, _, _, _ = MulpyRep.calcChopTimes(50.0, [50.0, 400.0], instpars, [17000], phaseoffset)
+        self.assertAlmostEqual(Eis_offset[5], 2.48991457, places=7)
+        self.assertAlmostEqual(Eis_offset[6], 2.10994937, places=7)
+        self.assertAlmostEqual(Eis_offset[7], 1.81075983, places=7)
+
 
 if __name__ == "__main__":
     unittest.main()

scripts/pychop/Instruments.py

@@ -135,6 +135,7 @@ class ChopperSystem(object):
         "flux_ref_slot",
         "flux_ref_freq",
         "frequency_names",
+        "phaseOffset",
     ]
 
     def __init__(self, inval=None):
@@ -168,6 +169,7 @@ def _parse_choppers(self):
         self.isFermi = False
         self.isPhaseIndependent = []
         self.defaultPhase = []
+        self.phaseOffset = None
         self.phaseNames = []
         for idx, chopper in enumerate(self.choppers):
             self.distance.append(chopper["distance"])
@@ -190,6 +192,8 @@ def _parse_choppers(self):
                 self.numDisk.append(2 if ("isDouble" in chopper and chopper["isDouble"]) else 1)
                 self.isPhaseIndependent.append(True if ("isPhaseIndependent" in chopper and chopper["isPhaseIndependent"]) else False)
                 self.defaultPhase.append(chopper["defaultPhase"] if "defaultPhase" in chopper else 0)
+                if "phaseOffset" in chopper:
+                    self.phaseOffset = chopper["phaseOffset"]
                 self.phaseNames.append(chopper["phaseName"] if "phaseName" in chopper else "Chopper %d phase delay time" % (idx))
         if not any(self.slot_ang_pos):
             self.slot_ang_pos = None
@@ -415,7 +419,9 @@ def _MulpyRepDriver(self, Ei_in=None, calc_res=True):
         """Private method to calculate resolution for given Ei from chopper opening times"""
         Ei = _check_input(self, Ei_in)
         if "_saved_state" not in self.__dict__ or (self._saved_state[0] != self._get_state(Ei)):
-            Eis, all_times, chop_times, lastChopDist, lines = MulpyRep.calcChopTimes(Ei, self._long_frequency, self._instpar, self.phase)
+            Eis, all_times, chop_times, lastChopDist, lines = MulpyRep.calcChopTimes(
+                Ei, self._long_frequency, self._instpar, self.phase, self.phaseOffset
+            )
             Eis, lines = self._removeLowIntensityReps(Eis, lines, Ei)
             self._saved_state = [self._get_state(Ei), Eis, chop_times, lastChopDist, lines, all_times]
         else:
@@ -645,7 +651,7 @@ def getMeasuredFlux(self, Ei):
         if not hasattr(self, "flux_interp"):
             raise AttributeError("This instrument does not have a table of measured flux")
         wavelength = [min(max(l, self.fmn), self.fmx) for l in np.sqrt(E2L / np.array(Ei if hasattr(Ei, "__len__") else [Ei]))]
-        return self.flux_interp(wavelength)
+        return self.flux_interp(wavelength[0])
 
     @property
     def theta_m(self):
@@ -995,7 +1001,6 @@ def calculate(cls, *args, **kwargs):
             except TypeError:
                 etrans = np.asfarray(etrans)
         res = obj.getResolution(etrans)
-
         if return_polynomial:
 
             def cubic(x, x_0, x_1, x_2, x_3):

scripts/pychop/MulpyRep.py

@@ -280,7 +280,7 @@ def calcFlux(Ei, freq1, percent, slot):
     return flux
 
 
-def calcChopTimes(efocus, freq, instrumentpars, chop2Phase=5):
+def calcChopTimes(efocus, freq, instrumentpars, chop2Phase=5, phaseOffset=None):
     """
     A method to calculate the various possible incident energies with a given chopper setup on LET.
     The window of energy transfers plotted is 85% by default.
@@ -310,6 +310,10 @@ def calcChopTimes(efocus, freq, instrumentpars, chop2Phase=5):
     ph_ind = np.array([False] * len(dist))
     if not len(ph_ind_v) == len(dist) and ph_ind_v:
         ph_ind[ph_ind_v] = True
+        # For Merlin, subtract the experimental offset of 4500
+        # This value can be set as phaseOffset in the input file for Merlin
+        if phaseOffset is not None:
+            chop2Phase[0] -= phaseOffset
         chop2Phase = phase = chop2Phase if hasattr(chop2Phase, "__len__") else [chop2Phase]
         if len(chop2Phase) != len(dist):
             if len(chop2Phase) == len(ph_ind_v):
@@ -378,4 +382,5 @@ def calcChopTimes(efocus, freq, instrumentpars, chop2Phase=5):
                 lines_all.append(line)
     # ok, now we know the possible neutron velocities. we now need their energies
     Ei = calcEnergy(lines_all, (dist[-1] + chop_samp))
+
     return Ei, chop_times, [chop_times[0][0], chop_times[-1][0]], dist[-1] - dist[0], lines_all

scripts/pychop/merlin.yaml

@@ -20,8 +20,9 @@ chopper_system:
       radius: 250               # Disk radius
       isDouble: False           # Is this a double disk system?
       isPhaseIndependent: True  # Is this disk to be phased independently?
-      defaultPhase: 1500        # What is the default phase for this disk (either a time in microseconds
+      defaultPhase: 17000        # What is the default phase for this disk (either a time in microseconds
                                 #   or a slot index [as a string] for the desired rep to go through)
+      phaseOffset: 4500         # Offset introduced after disk chopper controller update
       phaseName: 'Disk chopper phase delay time'
     -
       name: MERLIN Fermi