Fix bug in normalisation for POLARIS total scattering reduction

[RichardWaiteSTFC]

Describe the bug
S(Q) should tend to 1 at high Q if absolute units normalisation and conversion to S(Q) done correctly - but it doesn't!

FYI conversion from differential cross-section to S(Q)-1 is done

mantid/scripts/Diffraction/isis_powder/polaris_routines/polaris_algs.py

Lines 107 to 114 in 83a7b97

	# convert diff cross section to S(Q) - 1
	material_builder = MaterialBuilder()
	material_builder.setFormula(sample_details.material_object.chemical_formula)
	sample = material_builder.setNumberDensity(1.0).build() # number density not used in calculating xs (barns/atom)
	sample_total_scatter_cross_section = sample.totalScatterXSection()
	sample_coh_scatter_cross_section = sample.cohScatterXSection()
	focused_ws = focused_ws - sample_total_scatter_cross_section / (4 * math.pi)
	focused_ws = focused_ws * 4 * math.pi / sample_coh_scatter_cross_section

Absolute units is done here

mantid/scripts/Diffraction/isis_powder/routines/focus.py

Line 312 in 83a7b97

def _perform_absolute_normalization(spline, ws):

To Reproduce
Run the following script with this file s_of_q_minus_one_merged.nxs.txt

ws_sqm1 = LoadNexus(r"C:/Users/xhg73778/Documents/POLARIS/s_of_q_minus_one.nxs")
ws_sqm1 = ws_sqm1[1:] # for some reason group ws itself is first entry
ws_sqm1_merged = MatchAndMergeWorkspaces(InputWorkspaces=ws_sqm1, XMin=len(ws_sqm1)*[0.5], XMax=len(ws_sqm1)*[40], 
                                         CalculateScale=False)
                                         
ws_gr = PDFFourierTransform(InputWorkspace=ws_sqm1_merged_new, SofQType='S(Q)-1', PDFType='g(r)', 
                           DeltaR=None, RMax=20, Filter=False)  
                           
ws_sqm1_merged_bgsub = ws_sqm1_merged - ws_sqm1_merged.readY(0)[-1]
ws_gr_bgsub = PDFFourierTransform(InputWorkspace=ws_sqm1_merged_bgsub, SofQType='S(Q)-1', PDFType='g(r)', 
                                  DeltaR=None, RMax=20, Filter=False)  
                          

fig, axes = plt.subplots(dpi=156.0, edgecolor='#ffffff', num='ws_gr-1', subplot_kw={'projection': 'mantid'})
axes.plot(ws_gr, color='#1f77b4', label='ws_gr: spec 1', wkspIndex=0)
axes.plot(ws_gr_bgsub, color='#ff7f0e', label='ws_gr_bgsub: spec 1', wkspIndex=0)
axes.tick_params(axis='x', which='major', **{'gridOn': False, 'tick1On': True, 'tick2On': False, 'label1On': True, 'label2On': False, 'size': 6, 'tickdir': 'out', 'width': 1})
axes.tick_params(axis='y', which='major', **{'gridOn': False, 'tick1On': True, 'tick2On': False, 'label1On': True, 'label2On': False, 'size': 6, 'tickdir': 'out', 'width': 1})
axes.set_title('ws_gr')
axes.set_xlabel('Atomic Distance ($\\AA$)')
axes.set_ylabel('PDF')
axes.set_xlim([0.17274, 12.307])
axes.set_ylim([-2721.1, 2971.4])
legend = axes.legend(fontsize=8.0).set_draggable(True).legend
fig.show()

To run the POLARIS reduction follow instructions on this PR #39973

Expected behavior
S(Q) - 1 would tend to 0 at high-Q

[RichardWaiteSTFC]

Details of various normalisations can be found here in Pete Peterson's paper
https://onlinelibrary.wiley.com/doi/epdf/10.1107/S1600576720015630

In particular the conversion to S(Q) which is covered in Appendix C of the above paper (see Eq.41) and Appendix D (see Eq. 56 & 57)

[RichardWaiteSTFC]

FYI I tried this

@@ -108,12 +108,16 @@ def generate_ts_pdf(
     material_builder = MaterialBuilder()
     material_builder.setFormula(sample_details.material_object.chemical_formula)
     sample = material_builder.setNumberDensity(1.0).build()  # number density not used in calculating xs (barns/atom)
     sample_total_scatter_cross_section = sample.totalScatterXSection()
     sample_coh_scatter_cross_section = sample.cohScatterXSection()
-    focused_ws = focused_ws - sample_total_scatter_cross_section / (4 * math.pi)
-    focused_ws = focused_ws * 4 * math.pi / sample_coh_scatter_cross_section
+    # focused_ws = focused_ws - sample_total_scatter_cross_section / (4 * math.pi)
+    # focused_ws = focused_ws * 4 * math.pi / sample_coh_scatter_cross_section
+    bcoh_avg_sq = sample_coh_scatter_cross_section / (4 * np.pi)
+    btot_sq_avg = sample_total_scatter_cross_section / (4 * np.pi)
+    focused_ws = focused_ws - btot_sq_avg - bcoh_avg_sq
+    focused_ws = focused_ws * (1 / bcoh_avg_sq) - 1

Which I think is the same as Eq.41 in the Peterson paper (Assuming the material builder is doing what we think it is) - and gives similar results

It may be worth checking what MaterialBuilder is doing - it should be calculating quantities as in Table 3 of the errata to the Peterson paper published by Peterson and Keen
https://journals.iucr.org/j/issues/2021/05/00/gj5273/gj5273.pdf
i.e.
$<b_{tot}^2> = \sum_{i} c_i <b_i^2>$
$<b_{coh}^2> = \sum_{i} c_i <b_i^2>$
$<b_{coh}>^2 = (\sum_{i} c_i <b_i>)^2$
where
$c_i = N_i/N$ where $N_i$ is the number of atom species $N_i$ and $N$ is total number of atoms

[RichardWaiteSTFC]

FYI there is a mantid total scattering repo for Oak Ridge instruments (who are happy with their reduction)
This is the code for 'normalization' - this doesn't depend on the vanadium, but from the code it also looks like the workspace is already S(Q) at this point - so not sure at which point this is done.
https://github.com/neutrons/mantid_total_scattering/blob/7cff2bec01184f4bbf4d10c1bf8fd1665978bf87/total_scattering/reduction/normalizations.py#L47-L49
I think we need to go through this huge reduction function to see how we should be doing it!
https://github.com/neutrons/mantid_total_scattering/blob/7cff2bec01184f4bbf4d10c1bf8fd1665978bf87/total_scattering/reduction/total_scattering_reduction.py#L473

[RichardWaiteSTFC]

Script and config files etc. for BaTiO3 data

BTO_for_adrian.zip

[RichardWaiteSTFC]

Script and config files etc. for BaTiO3 data

BTO_for_adrian.zip

Also need script to create vanadium normalisation (you might need to change run-number to one of BTO runs)

van_spline_create.py.txt