From 28a82d8db73940fd0b8aef432c2edcc8edd3da49 Mon Sep 17 00:00:00 2001 From: Johannes Kasimir Date: Wed, 23 Jul 2025 17:27:48 +0200 Subject: [PATCH 01/10] feat: alternative beam center finder --- src/ess/sans/beam_center_finder.py | 54 ++++++++++++++++++++++++++++++ 1 file changed, 54 insertions(+) diff --git a/src/ess/sans/beam_center_finder.py b/src/ess/sans/beam_center_finder.py index ef9e82ef..b0ec3f67 100644 --- a/src/ess/sans/beam_center_finder.py +++ b/src/ess/sans/beam_center_finder.py @@ -37,6 +37,60 @@ def _xy_extrema(pos: sc.Variable) -> sc.Variable: return sc.concat([x_min, x_max, y_min, y_max], dim='extremes') +def _find_beam_center( + data, + sample_holder_radius=None, + sample_holder_arm_width=None, +): + if sample_holder_radius is None: + sample_holder_radius = sc.scalar(0.05, unit='m') + if sample_holder_arm_width is None: + sample_holder_arm_width = sc.scalar(0.015, unit='m') + + m = data.copy() + m.masks.clear() + s = m.bins.sum() + + for i in range(20): + c = (s.coords['position'] * sc.values(s)).sum() / sc.values(s).sum() + d = s.coords['position'] - c.data + + outer = 0.9 * min( + sc.abs(d.fields.x.min()), + sc.abs(d.fields.x.max()), + sc.abs(d.fields.y.min()), + sc.abs(d.fields.y.max()), + ) + s.masks['_outer'] = d.fields.x**2 + d.fields.y**2 > outer**2 + s.masks['_inner'] = d.fields.x**2 + d.fields.y**2 < sample_holder_radius**2 + + if i > 10: + s.coords['th'] = sc.where( + d.fields.x > sc.scalar(0.0, unit='m'), + sc.atan2(y=d.fields.y, x=d.fields.x), + sc.scalar(sc.constants.pi.value, unit='rad') + - sc.atan2(y=d.fields.y, x=-d.fields.x), + ) + h = s.drop_masks(['_arm'] if '_arm' in s.masks else []).hist(th=100) + th = s.coords['th'][np.argmin(h.values)] + + slope = sc.tan(th) / 2 + s.masks['_arm'] = ( + d.fields.y < slope * d.fields.x + sample_holder_arm_width + ) & (d.fields.y > slope * d.fields.x - sample_holder_arm_width) + return c.data + + +def beam_center_from_center_of_mass_alternative(workflow) -> BeamCenter: + try: + beam_center = workflow.compute(BeamCenter) + except sciline.UnsatisfiedRequirement: + beam_center = sc.vector([0.0, 0.0, 0.0], unit='m') + workflow[BeamCenter] = beam_center + data = workflow.compute(MaskedData[SampleRun]) + return _find_beam_center(data) + + def beam_center_from_center_of_mass(workflow: sciline.Pipeline) -> BeamCenter: """ Estimate the beam center via the center-of-mass of the data counts. From 1fe8ee0b186ce6d2895f085d85b414b835c4dcb5 Mon Sep 17 00:00:00 2001 From: Johannes Kasimir Date: Thu, 24 Jul 2025 09:23:00 +0200 Subject: [PATCH 02/10] fix: move arguments to public function + fix error in slope --- src/ess/sans/beam_center_finder.py | 26 +++++++++++++++----------- 1 file changed, 15 insertions(+), 11 deletions(-) diff --git a/src/ess/sans/beam_center_finder.py b/src/ess/sans/beam_center_finder.py index b0ec3f67..95fe74c0 100644 --- a/src/ess/sans/beam_center_finder.py +++ b/src/ess/sans/beam_center_finder.py @@ -39,14 +39,9 @@ def _xy_extrema(pos: sc.Variable) -> sc.Variable: def _find_beam_center( data, - sample_holder_radius=None, - sample_holder_arm_width=None, + sample_holder_radius, + sample_holder_arm_width, ): - if sample_holder_radius is None: - sample_holder_radius = sc.scalar(0.05, unit='m') - if sample_holder_arm_width is None: - sample_holder_arm_width = sc.scalar(0.015, unit='m') - m = data.copy() m.masks.clear() s = m.bins.sum() @@ -72,23 +67,32 @@ def _find_beam_center( - sc.atan2(y=d.fields.y, x=-d.fields.x), ) h = s.drop_masks(['_arm'] if '_arm' in s.masks else []).hist(th=100) - th = s.coords['th'][np.argmin(h.values)] + th = h.coords['th'][np.argmin(h.values)] - slope = sc.tan(th) / 2 + slope = sc.tan(th) s.masks['_arm'] = ( d.fields.y < slope * d.fields.x + sample_holder_arm_width ) & (d.fields.y > slope * d.fields.x - sample_holder_arm_width) return c.data -def beam_center_from_center_of_mass_alternative(workflow) -> BeamCenter: +def beam_center_from_center_of_mass_alternative( + workflow, + sample_holder_radius=None, + sample_holder_arm_width=None, +) -> BeamCenter: + if sample_holder_radius is None: + sample_holder_radius = sc.scalar(0.05, unit='m') + if sample_holder_arm_width is None: + sample_holder_arm_width = sc.scalar(0.02, unit='m') + try: beam_center = workflow.compute(BeamCenter) except sciline.UnsatisfiedRequirement: beam_center = sc.vector([0.0, 0.0, 0.0], unit='m') workflow[BeamCenter] = beam_center data = workflow.compute(MaskedData[SampleRun]) - return _find_beam_center(data) + return _find_beam_center(data, sample_holder_radius, sample_holder_arm_width) def beam_center_from_center_of_mass(workflow: sciline.Pipeline) -> BeamCenter: From 65ec69ad0d1d1f40b2cd3706673dbe17218054b2 Mon Sep 17 00:00:00 2001 From: Johannes Kasimir Date: Thu, 24 Jul 2025 15:56:33 +0200 Subject: [PATCH 03/10] docs: add docstrings explaining the procedure --- src/ess/sans/beam_center_finder.py | 64 ++++++++++++++++++++++++++++++ 1 file changed, 64 insertions(+) diff --git a/src/ess/sans/beam_center_finder.py b/src/ess/sans/beam_center_finder.py index 95fe74c0..98f387b8 100644 --- a/src/ess/sans/beam_center_finder.py +++ b/src/ess/sans/beam_center_finder.py @@ -42,6 +42,26 @@ def _find_beam_center( sample_holder_radius, sample_holder_arm_width, ): + ''' + Each iteration the center of mass of the remaining intensity is computed + and assigned to be the current beam center guess ``c``. + Then three symmetrical masks are created to make sure that the remaining intensity + distribution does not extend outside of the detector and that the sample holder + does not make the remaining intensity asymmetrical. + + The three masks are: + + - one "outer" circular mask with radius less than the minimal distance + from the current beam center guess to the border of the detector + - one "inner" circular mask with radius larger than the sample holder + - one "arm" rectangular mask with width wider than the sample holder arm + + The "outer" mask radius is found from the detector size. + The "inner" mask radius is supplied by the caller. + The "arm" mask slope is determined by the direction of minimum intensity + around the current beam center guess, the "arm" mask width is an argument + supplied by the user. + ''' m = data.copy() m.masks.clear() s = m.bins.sum() @@ -81,6 +101,50 @@ def beam_center_from_center_of_mass_alternative( sample_holder_radius=None, sample_holder_arm_width=None, ) -> BeamCenter: + """ + Estimate the beam center via the center-of-mass of the data counts. + + We are assuming the intensity distribution is symmetric around the beam center. + Even if the intensity distribution is symmetric around the beam center + the intensity distribution in the detector might not be, because + + - the detector has a finite extent, + - and there is a sample holder covering part of the detector. + + To deal with the limited size of the detector a mask can be applied that is small + enough so that the the remaining intensity is entirely inside the detector. + To deal with the sample holder we can mask the region of the detector that the + sample holder covers. + + But to preserve the symmetry of the intensity around the beam center the masks + also need to be symmetical around the beam center. + The problem is, the beam center is unknown. + However, if the beam center was known to us, and we applied symmetrical masks + that covered the regions of the detector where the intensity distribution is + asymmetrical, + then the center of mass of the remaining intensity would equal the beam center. + Conversely, if we apply symmetrical masks around a point that is not the beam center + the center of mass of the remaining intensity will (likely) not equal the original + point. + This suggests the beam center can be found using a fixed point iteration where each + iteration we + + 1. Compute the center of mass of the remaining intensity and assign it to be our + current estimate of the beam center. + 2. Create symmetrical masks around the current estimate of the beam center. + 3. Repeat from 1. until convergence. + + Parameters + ---------- + workflow: + The reduction workflow to compute MaskedData[SampleRun]. + + Returns + ------- + : + The beam center position as a vector. + """ + if sample_holder_radius is None: sample_holder_radius = sc.scalar(0.05, unit='m') if sample_holder_arm_width is None: From d36177cb35ff27cc3a30b9fc201767696327416c Mon Sep 17 00:00:00 2001 From: Johannes Kasimir Date: Thu, 24 Jul 2025 15:57:41 +0200 Subject: [PATCH 04/10] docs: fix --- src/ess/sans/beam_center_finder.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/ess/sans/beam_center_finder.py b/src/ess/sans/beam_center_finder.py index 98f387b8..4613be2e 100644 --- a/src/ess/sans/beam_center_finder.py +++ b/src/ess/sans/beam_center_finder.py @@ -60,7 +60,7 @@ def _find_beam_center( The "inner" mask radius is supplied by the caller. The "arm" mask slope is determined by the direction of minimum intensity around the current beam center guess, the "arm" mask width is an argument - supplied by the user. + supplied by the caller. ''' m = data.copy() m.masks.clear() From 37f6dc1220f0a5465e1391ebd6878d867b4322cd Mon Sep 17 00:00:00 2001 From: Johannes Kasimir Date: Wed, 12 Nov 2025 16:03:37 +0100 Subject: [PATCH 05/10] fix: add old beam center --- src/ess/sans/beam_center_finder.py | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/src/ess/sans/beam_center_finder.py b/src/ess/sans/beam_center_finder.py index 4613be2e..a5d3438c 100644 --- a/src/ess/sans/beam_center_finder.py +++ b/src/ess/sans/beam_center_finder.py @@ -93,6 +93,8 @@ def _find_beam_center( s.masks['_arm'] = ( d.fields.y < slope * d.fields.x + sample_holder_arm_width ) & (d.fields.y > slope * d.fields.x - sample_holder_arm_width) + + c.data.fields.z = sc.scalar(0.0, unit=c.data.unit) return c.data @@ -156,7 +158,10 @@ def beam_center_from_center_of_mass_alternative( beam_center = sc.vector([0.0, 0.0, 0.0], unit='m') workflow[BeamCenter] = beam_center data = workflow.compute(MaskedData[SampleRun]) - return _find_beam_center(data, sample_holder_radius, sample_holder_arm_width) + return ( + _find_beam_center(data, sample_holder_radius, sample_holder_arm_width) + + beam_center + ) def beam_center_from_center_of_mass(workflow: sciline.Pipeline) -> BeamCenter: From d889d911c587c237693013d625a3f40c6f704982 Mon Sep 17 00:00:00 2001 From: Johannes Kasimir Date: Wed, 12 Nov 2025 16:11:58 +0100 Subject: [PATCH 06/10] refactor: new domain types --- src/ess/sans/beam_center_finder.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/ess/sans/beam_center_finder.py b/src/ess/sans/beam_center_finder.py index dee354ef..2470412c 100644 --- a/src/ess/sans/beam_center_finder.py +++ b/src/ess/sans/beam_center_finder.py @@ -140,7 +140,7 @@ def beam_center_from_center_of_mass_alternative( Parameters ---------- workflow: - The reduction workflow to compute MaskedData[SampleRun]. + The reduction workflow to compute CorrectedDetector[SampleRun, Numerator]. Returns ------- @@ -158,7 +158,7 @@ def beam_center_from_center_of_mass_alternative( except sciline.UnsatisfiedRequirement: beam_center = sc.vector([0.0, 0.0, 0.0], unit='m') workflow[BeamCenter] = beam_center - data = workflow.compute(MaskedData[SampleRun]) + data = workflow.compute(CorrectedDetector[SampleRun, Numerator]) return ( _find_beam_center(data, sample_holder_radius, sample_holder_arm_width) + beam_center From 70cf3cf67a1fde5bd7d9581ac1fd5ddc0ff661ff Mon Sep 17 00:00:00 2001 From: Johannes Kasimir Date: Wed, 12 Nov 2025 16:17:44 +0100 Subject: [PATCH 07/10] docs: rename sample_holder to beam_stop --- src/ess/sans/beam_center_finder.py | 27 ++++++++++++--------------- 1 file changed, 12 insertions(+), 15 deletions(-) diff --git a/src/ess/sans/beam_center_finder.py b/src/ess/sans/beam_center_finder.py index 2470412c..12251f08 100644 --- a/src/ess/sans/beam_center_finder.py +++ b/src/ess/sans/beam_center_finder.py @@ -40,8 +40,8 @@ def _xy_extrema(pos: sc.Variable) -> sc.Variable: def _find_beam_center( data, - sample_holder_radius, - sample_holder_arm_width, + beam_stop_radius, + beam_stop_arm_width, ): ''' Each iteration the center of mass of the remaining intensity is computed @@ -78,7 +78,7 @@ def _find_beam_center( sc.abs(d.fields.y.max()), ) s.masks['_outer'] = d.fields.x**2 + d.fields.y**2 > outer**2 - s.masks['_inner'] = d.fields.x**2 + d.fields.y**2 < sample_holder_radius**2 + s.masks['_inner'] = d.fields.x**2 + d.fields.y**2 < beam_stop_radius**2 if i > 10: s.coords['th'] = sc.where( @@ -92,8 +92,8 @@ def _find_beam_center( slope = sc.tan(th) s.masks['_arm'] = ( - d.fields.y < slope * d.fields.x + sample_holder_arm_width - ) & (d.fields.y > slope * d.fields.x - sample_holder_arm_width) + d.fields.y < slope * d.fields.x + beam_stop_arm_width + ) & (d.fields.y > slope * d.fields.x - beam_stop_arm_width) c.data.fields.z = sc.scalar(0.0, unit=c.data.unit) return c.data @@ -101,8 +101,8 @@ def _find_beam_center( def beam_center_from_center_of_mass_alternative( workflow, - sample_holder_radius=None, - sample_holder_arm_width=None, + beam_stop_radius=None, + beam_stop_arm_width=None, ) -> BeamCenter: """ Estimate the beam center via the center-of-mass of the data counts. @@ -148,10 +148,10 @@ def beam_center_from_center_of_mass_alternative( The beam center position as a vector. """ - if sample_holder_radius is None: - sample_holder_radius = sc.scalar(0.05, unit='m') - if sample_holder_arm_width is None: - sample_holder_arm_width = sc.scalar(0.02, unit='m') + if beam_stop_radius is None: + beam_stop_radius = sc.scalar(0.05, unit='m') + if beam_stop_arm_width is None: + beam_stop_arm_width = sc.scalar(0.02, unit='m') try: beam_center = workflow.compute(BeamCenter) @@ -159,10 +159,7 @@ def beam_center_from_center_of_mass_alternative( beam_center = sc.vector([0.0, 0.0, 0.0], unit='m') workflow[BeamCenter] = beam_center data = workflow.compute(CorrectedDetector[SampleRun, Numerator]) - return ( - _find_beam_center(data, sample_holder_radius, sample_holder_arm_width) - + beam_center - ) + return _find_beam_center(data, beam_stop_radius, beam_stop_arm_width) + beam_center def beam_center_from_center_of_mass(workflow: sciline.Pipeline) -> BeamCenter: From 3cb087fbce74d9116af080416f71f4e4a7048040 Mon Sep 17 00:00:00 2001 From: Johannes Kasimir Date: Mon, 15 Dec 2025 10:04:40 +0100 Subject: [PATCH 08/10] refactor: variable naming --- src/ess/sans/beam_center_finder.py | 77 ++++++++++++++++++++---------- 1 file changed, 51 insertions(+), 26 deletions(-) diff --git a/src/ess/sans/beam_center_finder.py b/src/ess/sans/beam_center_finder.py index 12251f08..24f415ea 100644 --- a/src/ess/sans/beam_center_finder.py +++ b/src/ess/sans/beam_center_finder.py @@ -45,7 +45,7 @@ def _find_beam_center( ): ''' Each iteration the center of mass of the remaining intensity is computed - and assigned to be the current beam center guess ``c``. + and assigned to be the current beam center guess. Then three symmetrical masks are created to make sure that the remaining intensity distribution does not extend outside of the detector and that the sample holder does not make the remaining intensity asymmetrical. @@ -63,40 +63,65 @@ def _find_beam_center( around the current beam center guess, the "arm" mask width is an argument supplied by the caller. ''' - m = data.copy() - m.masks.clear() - s = m.bins.sum() + events = data.copy() + events.masks.clear() + intensity = events.bins.sum() for i in range(20): - c = (s.coords['position'] * sc.values(s)).sum() / sc.values(s).sum() - d = s.coords['position'] - c.data - + # Average position, weighted by intensity. + center = ( + intensity.coords['position'] * sc.values(intensity) + ).sum() / sc.values(intensity).sum() + distance_to_center = intensity.coords['position'] - center.data + + # Outer radius of annulus around center. + # Defined so that the annulus does not go outside of the bounds of the detector. outer = 0.9 * min( - sc.abs(d.fields.x.min()), - sc.abs(d.fields.x.max()), - sc.abs(d.fields.y.min()), - sc.abs(d.fields.y.max()), + sc.abs(distance_to_center.fields.x.min()), + sc.abs(distance_to_center.fields.x.max()), + sc.abs(distance_to_center.fields.y.min()), + sc.abs(distance_to_center.fields.y.max()), + ) + intensity.masks['outer'] = ( + distance_to_center.fields.x**2 + distance_to_center.fields.y**2 > outer**2 + ) + # Inner radius defined by size of beam stop. + intensity.masks['inner'] = ( + distance_to_center.fields.x**2 + distance_to_center.fields.y**2 + < beam_stop_radius**2 ) - s.masks['_outer'] = d.fields.x**2 + d.fields.y**2 > outer**2 - s.masks['_inner'] = d.fields.x**2 + d.fields.y**2 < beam_stop_radius**2 + # Iterate without the arm mask a few times to settle near the center + # before introducing the arm mask. if i > 10: - s.coords['th'] = sc.where( - d.fields.x > sc.scalar(0.0, unit='m'), - sc.atan2(y=d.fields.y, x=d.fields.x), + # Angle between +x and distance_to_center. + intensity.coords['theta'] = sc.where( + distance_to_center.fields.x > sc.scalar(0.0, unit='m'), + sc.atan2(y=distance_to_center.fields.y, x=distance_to_center.fields.x), sc.scalar(sc.constants.pi.value, unit='rad') - - sc.atan2(y=d.fields.y, x=-d.fields.x), + - sc.atan2( + y=distance_to_center.fields.y, x=-distance_to_center.fields.x + ), + ) + intensity_over_angle = intensity.drop_masks( + ['arm'] if 'arm' in intensity.masks else [] + ).hist(theta=100) + # Assume angle of arm coincides with intensity minimum + arm_angle = intensity_over_angle.coords['theta'][ + np.argmin(intensity_over_angle.values) + ] + + slope = sc.tan(arm_angle) + intensity.masks['arm'] = ( + distance_to_center.fields.y + < slope * distance_to_center.fields.x + beam_stop_arm_width + ) & ( + distance_to_center.fields.y + > slope * distance_to_center.fields.x - beam_stop_arm_width ) - h = s.drop_masks(['_arm'] if '_arm' in s.masks else []).hist(th=100) - th = h.coords['th'][np.argmin(h.values)] - - slope = sc.tan(th) - s.masks['_arm'] = ( - d.fields.y < slope * d.fields.x + beam_stop_arm_width - ) & (d.fields.y > slope * d.fields.x - beam_stop_arm_width) - c.data.fields.z = sc.scalar(0.0, unit=c.data.unit) - return c.data + center.data.fields.z = sc.scalar(0.0, unit=center.data.unit) + return center.data def beam_center_from_center_of_mass_alternative( From c8b4fb76cfbc50f8247a4b01c705f98206a76b53 Mon Sep 17 00:00:00 2001 From: Johannes Kasimir Date: Mon, 15 Dec 2025 10:07:59 +0100 Subject: [PATCH 09/10] docs: replace sample holder with beam stop --- src/ess/sans/beam_center_finder.py | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/src/ess/sans/beam_center_finder.py b/src/ess/sans/beam_center_finder.py index 24f415ea..30578ba2 100644 --- a/src/ess/sans/beam_center_finder.py +++ b/src/ess/sans/beam_center_finder.py @@ -47,15 +47,15 @@ def _find_beam_center( Each iteration the center of mass of the remaining intensity is computed and assigned to be the current beam center guess. Then three symmetrical masks are created to make sure that the remaining intensity - distribution does not extend outside of the detector and that the sample holder + distribution does not extend outside of the detector and that the beam stop does not make the remaining intensity asymmetrical. The three masks are: - one "outer" circular mask with radius less than the minimal distance from the current beam center guess to the border of the detector - - one "inner" circular mask with radius larger than the sample holder - - one "arm" rectangular mask with width wider than the sample holder arm + - one "inner" circular mask with radius larger than the beam stop + - one "arm" rectangular mask with width wider than the beam stop arm The "outer" mask radius is found from the detector size. The "inner" mask radius is supplied by the caller. @@ -137,12 +137,12 @@ def beam_center_from_center_of_mass_alternative( the intensity distribution in the detector might not be, because - the detector has a finite extent, - - and there is a sample holder covering part of the detector. + - and there is a beam stop covering part of the detector. To deal with the limited size of the detector a mask can be applied that is small enough so that the the remaining intensity is entirely inside the detector. - To deal with the sample holder we can mask the region of the detector that the - sample holder covers. + To deal with the beam stop we can mask the region of the detector that the + beam stop covers. But to preserve the symmetry of the intensity around the beam center the masks also need to be symmetical around the beam center. @@ -483,8 +483,8 @@ def beam_center_from_iofq( This was, however, not sufficient in cases where masks are applied to the detector pixels. It is indeed very common to mask broken pixels, as well as the region of - the detector where the sample holder is casting a shadow. - Such a sample holder will not appear in all 4 quadrants, and because it spans a + the detector where the beam stop is casting a shadow. + Such a beam stop will not appear in all 4 quadrants, and because it spans a range of scattering (:math:`2{\\theta}`) angles, it spans a range of :math:`Q` bins. All this means that we in fact need to perform a reduction as close as possible to From a32b7ebec6fe1e41253a027c7298e3c052e05445 Mon Sep 17 00:00:00 2001 From: Johannes Kasimir Date: Mon, 15 Dec 2025 13:00:50 +0100 Subject: [PATCH 10/10] test: result is unchanged --- tests/loki/iofq_test.py | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/tests/loki/iofq_test.py b/tests/loki/iofq_test.py index 08f3375d..c40fdfa5 100644 --- a/tests/loki/iofq_test.py +++ b/tests/loki/iofq_test.py @@ -283,6 +283,18 @@ def test_beam_center_from_center_of_mass_is_close_to_verified_result(): assert sc.allclose(center, reference) +def test_beam_center_from_center_of_mass_alternative_is_close_to_verified_result(): + pipeline = make_workflow(no_masks=False) + pipeline = sans.with_pixel_mask_filenames( + pipeline, loki.data.loki_tutorial_mask_filenames() + ) + center = sans.beam_center_finder.beam_center_from_center_of_mass_alternative( + pipeline + ) + reference = sc.vector([-0.0248743, -0.0166965, 0], unit='m') + assert sc.allclose(center, reference) + + def test_phi_with_gravity(): pipeline = make_workflow() pipeline[BeamCenter] = _compute_beam_center()