search

scipp / scippneutron

Neutron scattering toolkit built using scipp for Data Reduction. Not facility or instrument specific.
https://scipp.github.io/scippneutron/
BSD 3-Clause "New" or "Revised" License
4 stars 3 forks source link

Two theta gravity #514

Closed jl-wynen closed 2 months ago

jl-wynen commented 3 months ago

Fixes https://github.com/scipp/esssans/issues/90

~This is a draft because it does not include the reflectometry version of this function yet. But please review so I know whether the general approach is good.~

Reflectometry code postponed to later work.

~This version does not expose a way of computing unit vectors of the 'beam aligned coord system'. This would be needed to compute phi without gravity in SANS. But I don't know how to even define the coord system in that case. And the current implementation in ESSsans is not entirely clean because it has functions in the graph that depend on gravity.~

The basis vectors are now available. Projecting scattered_beam onto them not currently.

nvaytet commented 3 months ago

And the current implementation in ESSsans is not entirely clean because it has functions in the graph that depend on gravity.

I think we need gravity because we need to know which way is pointing up/down. Whether we include a correction for gravity in the computation of two_theta is separate, it depends on whether the function that computes two_theta takes gravity as an input or not?

jl-wynen commented 3 months ago

I think we need gravity because we need to know which way is pointing up/down.

True.

Whether we include a correction for gravity in the computation of two_theta is separate, it depends on whether the function that computes two_theta takes gravity as an input or not?

Technically, yes. This differs a little from other cases, though. Consider your proposal for ESSdiffraction where the presence of a TwoThetaBins param triggers focussing into two theta bins. So it seems that the presence of a gravity parameter would trigger gravity correction.

But we also need gravity to define a coordinate system. I currently think that it would be better to explicitly provide the basis vectors or at least an up vector instead of gravity.

jl-wynen commented 3 months ago

Do you have any comments about performance, especially memory usage? Compared to the implementation in ESSsans.

nvaytet commented 3 months ago

Do you have any comments about performance, especially memory usage? Compared to the implementation in ESSsans.

Isn't it trying to do the same (or pretty close to it?). One thing is that in essans we still need to define cylindrical_x and cylindrical_y for the beam center finder. How would the changes here fit in with that?

jl-wynen commented 3 months ago

Isn't it trying to do the same (or pretty close to it?).

Close. It differs in the definition for backscattering. The implementation here matches the definition in the regular two_theta function. This is why my implementation ultimately calls atan2 instead of asin.

One thing is that in essans we still need to define cylindrical_x and cylindrical_y for the beam center finder. How would the changes here fit in with that?

This could be done by adding a function / functions to compute those coordinates. They could then be inputs to scattering_angles_with_gravity. I did the computation inline here so that x can be dropped early to reduce memory usage.

Alternatively, we can provide a function to compute the basis vectors. (That is, make _beam_aligned_unit_vectors public.) The beam centre finder would then have to project scattered_beam onto those vectors itself and thus repeat that computation. But that is already the case in ESSsans.

Speaking of this, is the check in https://github.com/scipp/scippneutron/blob/e8319b850d4db4db9750cb12a18bd2fd7108faaf/src/scippneutron/conversion/beamline.py#L339-L347 too strict for the beam centre finder? Or for other applications? ESSsans simply assumed that gravity is orthogonal to the incident beam.

jl-wynen commented 3 months ago

@nvaytet This should be ready now. I updated the image and made beam_aligned_unit_vectors public. Should the computation of phi without gravity in ESSsans also be moved here?

nvaytet commented 3 months ago

is the check too strict for the beam centre finder? Or for other applications? ESSsans simply assumed that gravity is orthogonal to the incident beam.

I guess it depends when it would be enforced. For example, in relfectometry, the beam would not necessarily be horizontal for some instruments. For the beam center finder, you do get shifts that are of several millimeters, so the check could be too strict. We should try it.

nvaytet commented 3 months ago

Should the computation of phi without gravity in ESSsans also be moved here?

Could do, but leave for another PR?

nvaytet commented 3 months ago

I updated the image

Figures look good :+1:

jl-wynen commented 3 months ago

is the check too strict for the beam centre finder? Or for other applications? ESSsans simply assumed that gravity is orthogonal to the incident beam.

I guess it depends when it would be enforced. For example, in relfectometry, the beam would not necessarily be horizontal for some instruments. For the beam center finder, you do get shifts that are of several millimeters, so the check could be too strict. We should try it.

Ok. I will make a banch in ESSsans that uses the new code here to try it out.

jl-wynen commented 3 months ago

It works in ESSsans: https://github.com/scipp/esssans/pull/143

nvaytet commented 3 months ago

It works in ESSsans: scipp/esssans#143

By "works", did you mean that the old and new results are identical?

jl-wynen commented 3 months ago

The problem

There is indeed a discrepancy in the result. For example, take compute_Q of the SANS2d workflow. It produces CleanQ[SampleRun, Numerator]. Modifying compute_Q to convert to two_theta instead of Q and histogramming the result in two_theta, I get the following. The plot shows sc.values((new-old)/old) where new is the proposed code and old is the current main of ESSsans. two_theta_discrepancy

Looking at the two_theta coordinate in isolation, I get a max difference of ~1e-5 per spectrum-wavelength bin.

The cuprit

The difference stems from using atan2 in https://github.com/scipp/scippneutron/blob/33398a1f9d0280babfae3d33fc781d17d45fcb07/src/scippneutron/conversion/beamline.py#L516 where the implementation on ESSmain uses asin: https://github.com/scipp/esssans/blob/4d51289fe1205c0f7b113a52a944cccd61020ecd/src/ess/sans/conversions.py#L113 More concretely:

# proposed
two_theta = sc.atan2(y=sc.sqrt(x**2 + y**2), x=z)
# ESSsans main
two_theta = sc.asin(sc.sqrt(x**2 + y**2) / L2)

Using the asin-based code in the PR makes the SANS2D results (final and intermediate) agree perfectly with the implementation on ESSsans main.

I don't currently understand why this makes such a large difference. Naively, I would expect atan2 to have better precision than asin. But I can't tell whether this is the case here without further investigation.

jl-wynen commented 2 months ago

@nvaytet Can you take another look?