mantidproject / mantid

Main repository for Mantid code
https://www.mantidproject.org
GNU General Public License v3.0
210 stars 122 forks source link

Extend spherical absorption correction in AnvredCorrections to higher muR #31342

Closed RichardWaiteSTFC closed 2 years ago

RichardWaiteSTFC commented 3 years ago

AnvredCorrections corrects for the absorption of a spherical sample using a cubic fit to the A* values tabulated as a function of muR in Dwiggins (1975) [1] - for 5 degrees steps in theta. However, these values are limited to muR < 2.5 - however on WISH for some samples at long wavelength they exceed this limit.

In PR #31343 extrapolation of the cubic fit to Astar was allowed, however there exist values tabulated for muR>2.5 published in Weber (1967) [2] (in German!) which differs from the extrapolated Astar (although fractionally this is perhaps not too significant).

The below graph shows A* for theta=10 degrees. image

Although the values in [1] are slightly more accurate than [2] - the difference is typically much smaller the deviation from the cubic fit (the residuals seems to show a lot of structure so perhaps a higher order polynomial should be used?). One option is to replace the fit to [1] with a fit to [2] which goes up to muR<5. Another is just to use linear interpolation of [2] - no fit (interpolation is already required for different theta values).

[1] Dwiggins, C. W. "Rapid calculation of X-ray absorption correction factors for spheres to an accuracy of 0.05%." Acta Crystallographica Section A: Crystal Physics, Diffraction, Theoretical and General Crystallography 31.3 (1975): 395-396.

[2] Weber, K. "Eine neue Absorptionsfaktortafel fur kugelformige Proben." Acta Crystallographica Section B: Structural Crystallography and Crystal Chemistry 25.6 (1969): 1174-1178.

RichardWaiteSTFC commented 2 years ago

Found an error in the cubic polynomial coefficients for one theta value (but due to interpolation it would effect 45 < theta < 55)

image

RichardWaiteSTFC commented 2 years ago

The above is not enough to explain the discrepancy between the Monte-Carlo/Cubic mesh results and the AnvredCorrection result (see graph below - have checked both were converged). The output of AnvredCorrection appears to be correct - i.e. it agrees with Dwiggins [1]. I checked it against the parameterisation by Rouse et al. 1970 which gives transmission (1/A*) as a function of muR and theta - see Eq.2 in [3].

As pointed out in [3] - it is actually better to interpolate A* (in that paper they use Weber's values [2]) using sin(th/2)^2 - I can give this a go (note the results in [3] are only for muR<=1 - but perhaps including higher order terms in muR I can fit a greater range) [3] https://doi.org/10.1107/S0567739470001687 image

CreateSampleWorkspace(OutputWorkspace='ws', Function='Powder Diffraction', 
    XMin=2000, XMax=20000, BinWidth=20, NumBanks=1, BankPixelWidth=1, PixelSpacing=0.02,
    WorkspaceType="Histogram")
EditInstrumentGeometry(Workspace='ws', L2=0.5, Polar=100)
ConvertUnits(InputWorkspace='ws', OutputWorkspace='ws', Target='Wavelength')
shape = '''<sphere id="V-sphere">
    <centre x="0.0"  y="0.0" z="0.0" />
    <radius val="0.0025"/>
    </sphere>'''
SetSample('ws', Geometry={'Shape': 'CSG', 'Value': shape},
          Material={'ChemicalFormula': 'V0.95-Nb0.05', 'SampleNumberDensity': 0.0119})
SphericalAbsorption(InputWorkspace='ws', OutputWorkspace='attenuation_SphericalAbsorption')
w = MonteCarloAbsorption(InputWorkspace='ws', OutputWorkspace='attenuation_MonteCarlo', EventsPerPoint=1000)

wl = w.dataX(0)[:-1] + np.diff(w.readX(0)[:2])
mat = w.sample().getMaterial()
R = 100*(w.sample().getShape().volume()/(4*np.pi/3))**(1/3)
muR = np.array([(mat.totalScatterXSection() + mat.absorbXSection(wl_i))*mat.numberDensity*R for wl_i in wl]) # 
tth = w.spectrumInfo().twoTheta(0)
a1 = 1.5108
b1 = -0.1315
a2 = -0.0951
b2 = -0.2898
T = np.exp(-(a1 + b1*(np.sin(tth/2)**2))*muR - (a2 + b2*(np.sin(tth/2)**2))*(muR**2));
fig, ax = plt.subplots( subplot_kw={'projection': 'mantid'})
ax.plot(wl, T, '.', label='Rouse et al (1970)')
ax.set_xlabel('Wavelength')
ax.set_ylabel('Transmission')
ax.plot(mtd['attenuation_SphericalAbsorption'], distribution=True)
ax.plot(w, distribution=True)
ax.legend()
fig.show()
RichardWaiteSTFC commented 2 years ago

The agreement between Anvred and MonteCarlo seems pretty resonable ~5% even at larger muR - here I have increased the number density of a 2.5mm diameter sphere of Vanadium to change the muR range (note Rouse is only for muR <= 1) image