PowderN remove erroneous order parameter - #2622
Conversation
…e with the total cross section, not only the absorption cross section
…erefore it should be the total scattering cross section of the powder, and not only the incoherent xsect, no matter the order
|
Screenshot of relevant part of PowderN mcdoc page: 
Single Crystals documentation was not changed. |
|
@willend I am trying to run File "mcstas-dev/share/mcstas/tools/Python/mctest/mctest.py", line 979, in Have you seen this error before / know the work around? (The same command works for Single_crystal) |
|
The PowderN bug makes sense. Why was there an order parameter in PowderN, thought it was only single scattering? Was it done in preparation to adding multiple scattering to the component? I don't think it's wrong to only apply the absorption cross section for the exit path in Single_crystal when order is reached. After a scattering (which is known to have happened since the order is reached), the new wavevector will match the bragg condition of the crystal plane it already scattered in, and have a very high probability to scatter again, reducing the weight to almost 0 if it was used for the remaining path. What would happen if order was larger is the ray reflecting back and forth on that same reflection (and maybe more if it happens that more can be reached), so it would still have a big probability to continue in the direction order 1 predicted, but it would of course be attenuated by the path length, which is what the current order=1 approximates (the real path length would be slightly longer due to the zig-zag path). |
I also ran into this while trying to perform tests yesterday and did not figure it out. Will take a look soon. |
No. It was done, because I had misunderstood that attenuating by the total scattering cross section was a way to account for multiple scattering. (Which it isn't)
Hmmm. I don't quite think I am following you here. The way I see it, when order is reached, we make a monte carlo choice, to set the probability of transmitting to 1. The probability to actually transmit is exp(-total sigma * l), and therefore the factor we should multiply onto the weight is exactly exp(-total sigma * l) and not exp(-abs * l). With that said, my fix is not good enough, since the total scattering cross section that we access has not yet been modified for the direction of the neutron. |
Ah okay, so the order parameter that existed in PowderN had a different meaning than Single_crystal order.
The scattering back and forth on the same reflection is quite common and the current approximation doesn't quite follow it either, as there would only be a 50% chance of the ray leaving with that direction. I think your suggested change would be more accurate (order=1 being a good approximation to order=0) for very small crystals, while for larger crystals I would guess the current approximation would be closer. Could do some tests similar to PR #1767 |
|
@Lomholy @mads-bertelsen the mctest issue above seems associated with a bad state of the mctest on your branch... The one on main works nicely for me. I have therefore used the "Update branch" button on this PR to rebase against main. git pull, ./devel/bin/mccode-build-conda and you should hopefully be in a good place? Otherwise let me know and I can debug with you via zoom tomorrow. |
|
(as a result a new set of tests just started ...) |
Hi @willend I just tried to run the mctest with a git pull and conda build, as well as a deletion of the "old" mcstas stuff through
I will perform some tests like the linked PR, and we will see what happens! |
@Lomholy can you give me more details on version of OS, python, conda/micromamba etc? No chance to debug anything from what I currently have access to. |
Of Course. I am running a MacOs Tahoe 26.6.2, with a conda version 26.5.0 and a python version 3.14.4. Would a zoom call be easier? I have time until 11, and after 13. And I am doing it on an m4 max chip |
No Mads, It seems your original intuition was absolutely correct. I moved the attenuation down such that the beam is attenuated by the updated total cross section, and performed a test, as you can see at the bottom of this post, and it really does not work with the coherent cross section attenuating the beam. As you said, when we reflect the neutron is very likely to reflect again. This makes the importance choice of always leaving the crystal a quite poor approximation as any peak is simply reduced back in intensity. Therefore it is "not too bad" to simply attenuate with the absorption. However, is there then a reason not to attenuate with the incoherent cross section on the outgoing path? Below I made an example with an thin (0.1 mm) YBCO crystal with a fictitiously high incoherent cross section: 
Here the left hand image the total (outgoing) attenuation is applied, on the center figure both the incoherent and the absorption is applied, and on the right hand side only the absorption is applied. I am not quite sure if I have an argument for including the incoherent attenuation, other than the fact that it seems a better approximation, but neither do I have a proper mathematical argument for not including the total cross section (which provides a poor approximation). Here you have the same plot for a broader theta range: |
I encountered the bug too and thought it was just something with my build. Looked through it with help of AI and believe I have a good fix on PR #2624 |
Great work! My argument for not including incoherent scattering in the attenuation is that the neutrons coming from that reflection to the detector should be attenuated, but rays that are scattered elsewhere should also be incoherently scattered towards this peak, but since we don't simulate the added intensity from the incoherent elsewhere, we are canceling out the effect. You can add extra_order to have additional incoherent scattering after the base order is reached, then you get both the scattering and attention for very little performance cost. It can mean a lot in crystals with significant incoherent cross section. For some reason I can't download the graphs, but I do see including it seems to converge faster, so my intuition on not including incoherent seems wrong. Is it also better for a large crystal? |
@mads-bertelsen thanks. I think I'd rather do this one by hand. |
All good! Hope the PR helps identify the underlying issue |
mctest aborts a whole test run with
TypeError: %d format: a real number is required, not NoneType
in the "Display OK/FAILED" status line whenever it re-enters an
instrument test dir that already contains results from an earlier run.
Reported by users on PR #2622 (crash at the "Display FAILED" line).
How it happens:
- If the per-instrument test dir already exists - e.g. a re-run with
the same --testdir and --uid, which is exactly how the CI workflows
invoke mctest - shutil.copytree fails and mctest falls back to
InstrExampleTest.load() to rebuild the test objects from the JSON
files saved by the previous run.
- load() restored every stored field except "displaytime", leaving the
test object with displayed=<saved value> but displaytime=None. The
"binary already exists" branch of the compile loop skips the
mcdisplay step as well, so displaytime is never (re)measured.
- The run/status loop then formats test.displaytime with %d, which
raises TypeError and kills the run before any results are written.
The "Display FAILED" variant (the one seen on PR #2622) occurs when
the first run recorded displayed=false, e.g. because mcdisplay hit the
60 s --displaymax timeout - seen for MCPL instruments such as
ESS_BEER_MCPL / ESS_KVASIR / ESS_butterfly_MCPL_test.
Recreate the error (before this fix):
mctest --testdir tmp --instr ESS_KVASIR -n 1e4 --uid reproA
# 1st run: SUCCESS, saves results JSON
mctest --testdir tmp --instr ESS_KVASIR -n 1e4 --uid reproA
# 2nd run: "WARNING: Skipped ESS_KVASIR test - did ... exist
# already??" followed by the TypeError traceback above
The FAILED-branch variant is forced by making the display time out in
the first run, e.g. adding --displaymax 0.
Fix:
- InstrExampleTest.load() now restores displaytime from the saved JSON
(via obj.get, so JSON files written without the key still load).
- The Display OK/FAILED status lines fall back to 0 s when displaytime
is None, which also covers a binary present without saved display
time (e.g. a first run interrupted between compile and saving).
Verified: both re-run scenarios above now complete with SUCCESS, and
fresh mctest runs are unaffected.
Below I conducted the test again but changed the incoherent cross section a bit, and made the crystal thick. The convergence does not seem to be certain across parameters. Overall I think the approach with just setting the extra order high is probably better than this. I will remove the addition from this PR. 
Thank you both for swift responses and good discussions! |
|
The PowderN component lies within quite a lot of instruments, and I therefore only post screenshots of the failed tests. 
I will run it again at higher n count to see if the problem persists with the test_powders instrument. I am however uncertain as to how this change could have affected any display function, so the display errors are likely not from this PR. |
|
When running the Powders16 test again I get 88 % with ncount=1e7: 
And the display errors still apply. |
|
@Lomholy I will run a series of tests just with Test_Powders to to check if we need to (slightly) update the test value for _16... Am I correct that the PR title ended up being slightly off? (It looks like only PowderN was touched in the end?) |
Cool Thank you!
Yes, I will update it. Only PowderN was touched upon in the end :D |
Lomholy
changed the title
3 participants
Free-form text area
Please describe what your PR is adding in terms of features or bugfixes:
This PR fixes two minor bugs on the Single_crystal component, and the PowderN component.
The bug in the Single_Crystal component is, that when the number of scatterings is exceeded due to the order parameter, the ray is weighted with exp(-absl). This is wrong, since it should be the total cross section of events i.e exp(-totl).
The bug in the PowderN component is that the attenuation of the ray due to sampling the position is wrong. Currently when order is not set, then my_s is set to my_inc, and it should always be set to the total cross section.
I (embarrasingly enough) implemented this erroneous parameter myself, due to a misunderstanding of the internal workings of this component. I believed that my_s was set to the total cross section in order to mitigate multiple scattering, when it in fact was only there to correct for sampling uniformly instead of exponentially.
Development OS / boundary conditions
Please describe what OS you developed and tested your additions on, and if any special dependencies are required:
MacOS Tahoe 26.6.2
PR Checklist for contributing to McStas/McXtrace
For a coherent and useful contribution to McStas/McXtrace, please fill in relevant parts of the checklist:
My contribution includes patches to an existing component file
mcdocutility and rendered a reasonable documentation page for the component (please attach as screenshot in comments!)mctestutility to test one or more instruments making use of the component (please attachmcviewtestreport as screenshot in comments)mccode-clangformattool to apply the standard McCode component indentation schememcrun --c-lint"linter" and followed advice to remove most / all warnings that are raised