Update pychop definitions for SNS instruments

[mducle]

Description of work

When SNS instruments were added to PyChop around 2020 I used the original .yaml files rather than fine-tuned files which Jiao Lin had created during the work for this publication. So, the calculations in Mantid has actually been incorrect for 4 years!

During work for the resconvlib project @RastislavTuranyi found that Mantid calculations disagreed with the online calculator provided by SNS which is also meant to be using PyChop. Furthermore the Mantid calculations disagreed drastically with the measured data in the ARCS resolution paper.

The discrepancy is because the moderator parameters (amongst other things) were incorrect in the initial version added to Mantid.

Summary of work

This PR is to update these .yaml files to the optimised files in Jiao Lin's original repository, namely the files:

• ARCS/PyChop/arcs-09122018.yaml
• SEQUOIA/PyChop/sequoia-06212019.yaml

There is no associated issue.

Further detail of work

To test:

1. Run Mantid and open the PyChop interface (Direct->PyChop).
2. Select ARCS and type in an incident energy (say 100meV) and click calculate.
3. Go to the SNS resolution webpage; select ARCS and use the same parameters as in the Mantid/PyChop calculation, and check the calculations agree.
4. Repeat for Sequoia.

---

Reviewer

Please comment on the points listed below (full description).
Your comments will be used as part of the gatekeeper process, so please comment clearly on what you have checked during your review. If changes are made to the PR during the review process then your final comment will be the most important for gatekeepers. In this comment you should make it clear why any earlier review is still valid, or confirm that all requested changes have been addressed.

Code Review

• Is the code of an acceptable quality?
• Does the code conform to the coding standards?
• Are the unit tests small and test the class in isolation?
• If there is GUI work does it follow the GUI standards?
• If there are changes in the release notes then do they describe the changes appropriately?
• Do the release notes conform to the release notes guide?

Functional Tests

• Do changes function as described? Add comments below that describe the tests performed?
• Do the changes handle unexpected situations, e.g. bad input?
• Has the relevant (user and developer) documentation been added/updated?

Does everything look good? Mark the review as Approve. A member of @mantidproject/gatekeepers will take care of it.

Gatekeeper

If you need to request changes to a PR then please add a comment and set the review status to "Request changes". This will stop the PR from showing up in the list for other gatekeepers.

[SilkeSchomann]

@mducle Please add a release note in \docs\source\release\v6.12.0\Direct_Geometry\General\Bugfixes

[ajjackson]

I notice that in the .yaml for SEQUOIA, the comment on moderator parameters has not changed and indicates that B will be used directly (rather than as reciprocal).

I compared these comments a bit more closely with the code, and neither of them seems to be an accurate description of the tikeda function

mantid/scripts/pychop/Chop.py

Lines 116 to 129 in 6ee0b15

	def tikeda(S1, S2, B1, B2, Emod, Ei):
	"""
	! Calculates the moderator time width based on the Ikeda-Carpenter distribution
	"""
	Ei = np.array(Ei if hasattr(Ei, "__len__") else [Ei])
	sig = np.sqrt((S1 * S1) + ((S2 * S2 * 81.8048) / Ei))
	A = 4.37392e-4 * sig * np.sqrt(Ei)
	tausqr = []
	B = np.array([B1] * len(Ei))
	B[np.where(Ei > 130.0)] = B2
	R = np.exp(-Ei / Emod)
	tausqr = (3.0 / (A**2)) + (R * (2.0 - R)) / (B**2)
	# variance currently in mms**2. Convert to sec**2
	return (tausqr if len(tausqr) > 1 else tausqr[0]) * 1.0e-12

Specifically:

• The final expression is effectively $\tau^2 = 3 \tau_f^2 + R(2-R)\tau_s^2$; the components of $\tau$ are not reciprocals. (This makes sense as it means all the $\tau_i$ have the same units.) The new ARCS comment seems to be correct about the treatment of B1/B2 here.
• I can't figure out where the conversion 81.8048 comes from. It should be converting an energy in meV to a wavelength according to the Ikeda-Carpenter paper.
• The $R(2-R)$ group in the implementation looks a bit odd, and differs from the comments which describe $R/(2-R)$. Neither form obviously resembles the Ikeda-Carpenter function in the literature; I think there must be some algebraic simplification at work?

refs:

• Mantid docs
• Ikeda-Carpenter paper
• Paper on fitting Ikeda-Carpenter functions for instruments

[mducle]

81.81 is a "well-known" conversion factor for neutron conversion where E = 81.81 / lambda**2 when energy is in meV and wavelength is in Angstrom. It amounts to h/2m.

I'll see if I can dig out where the approximation comes from... it's probably more fully explained in the original Fortran code, which I have... somewhere...

Also the unit test is broken because it compares numerical values of the resolution and flux which is obviously now different (but correct!). I'll fix that and add release notes.

[AndreiSavici]

@ajjackson The 81.81 is the conversion from energy in meV to wavelength in Angstroms. $h^2/(2m_n\lambda^2)=\Delta E$ But this formula is in SI. We use $\lambda$ in Angstroms, and $\Delta E (meV)=\Delta E(J) 10^{-3}e$, where $e$ is the elementary charge. So the constant is

$$\frac{h^2}{2m_n e (10^{-10})^210^{-3}}$$

import scipy.constants as sc
lam2meV = sc.Planck**2/(2*sc.e*sc.m_n*1e-23)

[mducle]

@ajjackson The original Fortran code is in a gist here. Not sure it's super useful. I think that the value returned by the function, tausqr is variance of the distribution in eq 10' of Ikeda-Carpenter so it should be $\int \psi^2(v, t) dt$... I might see if I can plug that into Wolfram...

[mducle]

@AndreiSavici I've changed the .yaml files here a bit compared to the arcs-09122018.yaml and sequoia-06212019.yaml files in order to match what the website gives instead.

For ARCS I've included a scale factor of 600 for the flux and for SEQUOIA I've changed the look up table to include the detector efficiency and absorption corrections which the website code does.

So now Mantid-PyChop will give the same flux (in n/cm^2/s/MW) as the webpage.

@SilkeSchomann I've also changed the documentation to explicitly say that the flux values given for ISIS and SNS instruments are in different units. (Although in principle if both facilities are running at full power the flux should be in the same units of n/cm^2/s). I haven't changed the GUI - it still says n/cm^2/s - do you think that's needed for this?

[SilkeSchomann]

@AndreiSavici I've changed the .yaml files here a bit compared to the arcs-09122018.yaml and sequoia-06212019.yaml files in order to match what the website gives instead.

For ARCS I've included a scale factor of 600 for the flux and for SEQUOIA I've changed the look up table to include the detector efficiency and absorption corrections which the website code does.

So now Mantid-PyChop will give the same flux (in n/cm^2/s/MW) as the webpage.

@SilkeSchomann I've also changed the documentation to explicitly say that the flux values given for ISIS and SNS instruments are in different units. (Although in principle if both facilities are running at full power the flux should be in the same units of n/cm^2/s). I haven't changed the GUI - it still says n/cm^2/s - do you think that's needed for this?

@mducle I would change the GUI as well. A mismatch between GUI and documentation will probably confuse the user.