Skip to content

[BUG]: McStas: code written for 2.x in-place RESTORE_NEUTRON; an ABSORB after RESTORE is cancelled in 3.x #2939

Open

Description

@willend

Background

In McStas 2.x, RESTORE_NEUTRON(index, x, y, z, vx, vy, vz, t, sx, sy, sz, p) restored the state in place, immediately (mcrestore_neutron(mccomp_storein, index, &x, ...) in mcstas-r.h). Since 3.0 it only sets a flag (cogen.c.in: #define RESTORE_NEUTRON(_index, ...) _particle->_restore = _index;). The generated code then restores the component's entry state after TRACE and EXTEND, and particle_restore() also clears _absorbed:

class_X_trace(&_X_var, _particle);
if (_particle->_restore)
  particle_restore(_particle, &_particle_save);   // ... p->_absorbed=0; p->_restore=0;

Code written for the old semantics therefore behaves differently in 3.x:

  1. An ABSORB after a RESTORE_NEUTRON, in TRACE or in the instrument's EXTEND, is cancelled. The neutron continues downstream.
  2. Code that keeps going after RESTORE_NEUTRON works on the un-restored state, and everything it does (propagation, weight, velocity) is then discarded at the end of the component. SCATTERED stays set, though.
  3. The index argument is ignored. You can only go back to the current component's entry state, not to an earlier component.

Survey of origin/main (e2d7eb1), McStas side. The McXtrace side is a separate issue.

Reproducer

Restore_absorb_n.instr
DEFINE INSTRUMENT Restore_absorb_n(int mode=1)
TRACE
COMPONENT src = Source_simple(radius=0.01, dist=1, focus_xw=0.02, focus_yh=0.02, lambda0=4, dlambda=0.1)
AT (0,0,0) ABSOLUTE
COMPONENT mono = Monochromator_curved(width=0.01, height=0.01, NH=1, NV=1, DM=3.355)
WHEN (mode==1) AT (0.5,0,0.5) RELATIVE src ROTATED (0,45,0) RELATIVE src
EXTEND %{
  if (!SCATTERED) ABSORB;
%}
COMPONENT mnd = Monitor_nD(xwidth=0.001, yheight=0.001, options="x y")
WHEN (mode==2) AT (0.5,0,0.6) RELATIVE src
EXTEND %{
  ABSORB;
%}
COMPONENT mnda = Monitor_nD(xwidth=0.001, yheight=0.001, options="absorb x y")
WHEN (mode==3) AT (0.5,0,0.7) RELATIVE src
COMPONENT arm = Arm()
AT (0,0,0.8) RELATIVE src
EXTEND %{
  if (INSTRUMENT_GETPAR(mode)==4) ABSORB;
%}
COMPONENT det = PSD_monitor(xwidth=0.1, yheight=0.1, filename="det")
AT (0,0,1) RELATIVE src
END

Each case places a component off-beam and absorbs the neutrons it doesn't scatter, so det should get 0 everywhere. Results with 1e5 neutrons (mcstas 3.99.99, macOS arm64):

mode case events at det
1 Monochromator_curved missed, EXTEND if (!SCATTERED) ABSORB 49798 (all the restored ones)
2 Monitor_nD missed, EXTEND ABSORB 100000
3 Monitor_nD missed, options="absorb" 100000
4 control: Arm, EXTEND ABSORB 0

Affected: ABSORB after RESTORE in a component's TRACE

  • monitors/Monitor_nD.comp:666 and monitors/Monitor_nD_noacc.comp:659: RESTORE_NEUTRON(...); ABSORB; in the no-intersection branch when absorb is set. Because of the RESTORE, the absorb option has no effect on neutrons that miss (mode 3). The comment says the restore is only there "for correct mcdisplay".

Affected: code that relies on the restore happening in place

  • contrib/Multilayer_Sample.comp:127: in the incoherent branch, RESTORE_NEUTRON is followed by the whole V_sample-like scattering calculation and several ABSORBs. All of that is undone at the end of the component, and the ABSORBs are cancelled. Line 165 (the "leave alone" branch) is the same.
  • examples/TUDelft/SESANS_Delft/SESANS_Delft.instr:172: the trick Arm calls RESTORE_NEUTRON(INDEX_CURRENT_COMP-2, ...) in EXTEND, to send the pre-analyser state through ana2/det2. In 3.x this restores the Arm's own entry state, i.e. the neutron that has already passed ana at z = 4.5 m. Result: det2 is exactly 0 (1e6 neutrons; det = 2.26e-7). Only det is covered by %Example, so tests don't catch this.

Affected: instruments with EXTEND ABSORB after a component that RESTOREs

These components restore on a miss (or always, with restore_neutron=1), so for those neutrons the instrument's ABSORB is cancelled:

  • Monochromator_curved + if (!SCATTERED) ABSORB: ILL_IN4 Mono, ILL_IN13 IN13_deflector and Analyser, ILL_Lagrange Analyser, ILL_H5_new H521_SADAM_Monochromator, ILL_H521_NeXT NeXT_DCM_1, Test_focus mono_ideal, Test_TasReso mono and ana
  • Pol_mirror + if (!SCATTERED) ABSORB: SE_example and SE_example2, polariser2 and analyzer
  • Monitor_nD used as a beam dump (ABSORB; /* neutron ends on the ... sample */): misses survive and continue into the following branches. ILL_H5 (H53_IN14/IN16/D16/SADAM_Diffractometer, H511_IN15_Diffractometer, D22_Detector), ILL_H512_D22 D22_Detector, ILL_H5_new (H53_ThALES_Diffractometer, H53_ThALES_TransmitXY, H521_D16/SADAM_Diffractometer, H522_WASP_Diffractometer, H511_IN15_Diffractometer, D22_Detector)
  • Monitor_nD(restore_neutron=1) + EXTEND ABSORB: ILL_IN13 CaF2_transmitted (if (SCATTERED) ABSORB), ILL_IN5_Spots Det_sample_t (if (!SCATTERED) ABSORB)
  • Monitor_nD in a GROUP + if (!mono_index) ABSORB: TRIGA/RTP_SANS psd_reflect (needs checking together with the GROUP restore logic)

I have not quantified the impact per instrument. Whether the leaked neutrons reach a detector depends on the downstream geometry.

Not affected (as far as I can see)

  • RESTORE as the last action on a miss, with no EXTEND ABSORB in the instruments: Beamstop, PerfectCrystal, Pol_bender_tapering, NPI_tof_dhkl_detector, and so on.
  • Monitors that restore at the end of TRACE (restore_neutron=1), unless an instrument ABSORBs after them (see the list above).
  • if (SCATTERED) ABSORB after a component that doesn't restore on scatter, e.g. ILL_H521_NeXT H521_D16_Monochromator.

Possible directions (for discussion)

  • Per component, without touching the code generator: on a miss, undo the propagation by hand (save x, y, z, t before, and put them back) instead of RESTORE_NEUTRON, so an EXTEND ABSORB wins again. This is the approach taken for McXtrace Mirror_curved in the follow-up to Fix curved mirror #2936. Monitor_nD: drop the RESTORE before ABSORB. Multilayer_Sample: keep local copies instead of RESTORE. SESANS_Delft: rework the trick (store the state in USERVARS before ana, and reset it explicitly).
  • Centrally: skip particle_restore() when ABSORBED is set, so ABSORB always beats RESTORE. That is one line in the code generator and would fix all the EXTEND cases at once, but it is a code-generator change and changes numbers.

Either way this changes simulation results for the instruments listed, so %Example targets need re-checking.

Related: #1227 (RESTORE and USERVARS in McCode 3).

Survey and reproducers made with Claude Code.

Activity

  1. willend commented on Oct 30, 2026

    @willend
    ContributorAuthor

    McXtrace counterpart: #2940

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Metadata

Metadata

Assignees

No one assigned

    Labels

    No labels
    No labels

    Type

    No type

    Projects

    No projects

      Milestone

      No milestone

      Relationships

      None yet

      Development

      No branches or pull requests

      Issue actions