reflectometry / refl1d

1-D reflectometry fitting
https://refl1d.readthedocs.io/
Other
17 stars 24 forks source link

better support for Q R dR files #47

Open pkienzle opened 12 years ago

pkienzle commented 12 years ago

Many facilities produce Q, R, dR files with constant dQ/Q.

Include a function like read_QRdR with appropriate resolution options to make this easier.

For example

def read_QRdR(filename, dQoQ): Q,R,dR = numpy.loadtxt(filename).T probe = QProbe(Q, dQoQ*Q, data=(R,dR)) probe._sf_L = numpy.array([1.])

probe._sf_idx = numpy.zeros(len(Q),'int32')

probe.calc_Qo = probe.Q
return probe

pp = read_QRdR('barb1A.ref', dQoQ=0.01) mm = read_QRdR('barb1D.ref', dQoQ=0.01) probe = PolarizedNeutronQProbe([mm,None,None,pp])

pkienzle commented 4 years ago

load4 is now much more flexible, so you can indicate θ, Δθ, λ, Δλ directly in the file header (as angle, angular_resolution, wavelength, wavelength_resolution) or as parameters (T, dT, L, dL). These can be lists if they differ for each element, or even functions (dT=lambda T: ..., dL=lambda L: ...).

There is still no simple dQ for those datasets that mix both θ and λ, such as TOF at multiple angles or white beam reflectometers with wavelength, so the ticket needs to stay open.

pkienzle commented 3 years ago

In practice on Candor the resolution on the measured q points is tiny relative to the width of the q bins and the q resolution is approximately uniform rather than gaussian. Furthermore, given the spectrum we can use the average wavelength to estimate the average angle for the given q, and from that compute q offset as q' = q sin(<θ> + radians(δθ))/sin(<θ>) ≈ q(1 + δθ/<θ>). Getting extra fancy, θ offset for binned q values will broaden the q resolution, and this, too, can be approximated using the wavelength spectrum.

Similarly for sample broadening, but there we may need the nominal Δθ as well.

pkienzle commented 3 years ago

Demonstration of estimated vs actual q-offset using average of the wavelength spectrum:

λ = np.linspace(4, 6, 10)  # wavelength 4-6 Å
q = 0.2                # high q
δ = radians(0.1)  # 0.1° angle offset

# true q' after theta offset
π = np.pi
Eλ = mean(λ)
θ = np.arcsin(q*λ/(4*π))  # actual angle
qp = 4*π/λ*sin(θ+δ)
print(mean(qp), std(qp))
# 0.2044458970751779 0.0005822762660856065

# estimated q' after theta offset
Eθ = np.arcsin(q*Eλ/(4*π))  # angle estimated from average wavelength
print(q*sin(δ+Eθ)/sin(Eθ))
# 0.2043722730607561
print(q*(1+δ/Eθ)) # small angle approximation
# 0.20438185288153524
pkienzle commented 3 years ago

Using average wavelength for θ offset only works for a relatively narrow wavelength band.

Using λ ∈ [2.5, 17] with a Maxwell distribution peaked at 4 Å the approximation is not sufficient.

# Neutron spectrum from Maxwell distribution with 1/ν detection efficiency
# as given in Eq 9 of http://web.mit.edu/8.13/www/JLExperiments/JLExp38.pdf
def E(λ): return 81.82/λ**2
def ν(E): return 437*sqrt(E)
def rate(ν, νo): return ν**2*exp(-ν**2/νo**2)

λ = np.linspace(2.5, 17, 10)  # wavelength 2.5 - 17 Å
wλ = rate(ν(E(λ)), 1000)
q = 0.02
δ = radians(0.1)  # 0.1° angle offset

# true q' after theta offset
π = np.pi
Eλ = np.average(λ, weights=wλ)
θ = np.arcsin(q*λ/(4*π))  # actual angle
qp = 4*π/λ*sin(θ+δ)
Eqp = np.average(qp, weights=wλ)
Vqp = np.average((qp-Eqp)**2, weights=wλ)
print(Eqp, np.sqrt(Vqp))
# 0.024016859114951838 0.0022114219251348744

# estimated q' after theta offset
Eθ = np.arcsin(q*Eλ/(4*π))  # angle estimated from average wavelength
print(q*sin(δ+Eθ)/sin(Eθ))
# 0.023005158893326525
print(q*(1+δ/Eθ)) # small angle approximation
# 0.02300532602758869
bmaranville commented 3 years ago

The Maxwell distribution part may be close to what is seen in a real experiment, but I doubt that 1/nu describes the detection efficiency in the wavelength range we look at with real detectors, does it?

pkienzle commented 3 years ago

I don't know. Here's the "exercise for the student" following Eq 9:

Can you figure out the reason for using a “thin” counter rather than using one with higher efficiency of say 50%, aside from the fact that we don’t need the higher efficiency because the measured intensity is adequately high?

The issue is worse using ν³ instead of ν²

0.025234202764418703 0.002441826169924171
0.024016910178888518
0.024017043788838656

The real spectrum is going to be convoluted with proton pulse shape, and I'm sure many other refinements. We can try the calculation with a measured spectrum if you have one, but I don't think it will change the conclusion that we can't cheat with θ-offset when we have a broad spectrum with a strong peak.

bmaranville commented 3 years ago

agreed

pkienzle commented 3 years ago

Adjusting the correction slightly, the follow works very well for q ∈ [0.001, 0.3] and δθ ∈ [0.001º, 1°]:

Eλ = 1/average(1/λ, weights=wλ)
Vλ = 1/average((1/λ - 1/Eλ)**2, weights=wλ)
Eθ = arcsin(q*Eλ/(4*π))
qp = q*sin(np.radians(δ)+Eθ)/sin(Eθ)
Δqp_broadening = 4*π/sqrt(Vλ)*sin(δ)

with q' error less than 0.02% and Δq' error less than 0.7%.

The Δq broadening expression comes from the following:

q' = 4π/λ sin(θ+δθ) with sin θ = λq/4π for constant q
E[q'] = ⟨4π/λ sin(θ+δθ)⟩_λ
V[q'] = ⟨(4π/λ sin(θ+δθ) - E[q'])² ⟩_λ

sin(θ+δθ) ≈ sin θ + sin δθ
⇒ V[q'] ≈ ⟨(q + 4π/λ sin δθ - ⟨q + 4π/λ sin δθ⟩_λ)² ⟩_λ
        = (4π sin δθ)² ⟨(1/λ - ⟨1/λ⟩_λ)² ⟩_λ

From the plot below there is clearly a q dependent correction we could make for the small angle approximation but we are not computing the resolution accurately enough for it to matter.

Conclusion: given the mean and variance in the spectrum we can fit δθ offset even on stitched TOF data.

image

image

pkienzle commented 3 years ago

We still need an expression for sample broadening on stitched data. That is, if the sample is slightly convex and there is increased divergence, then how does this affect the ΔQ resolution?