Epithermal and fast neutron interactions / sapphire filter

[9klaud]

Hi there!

As part of my PhD project, I developed a validated Monte Carlo (Geant4), CAD-based model of a thermal neutron beamline and I plan to implement a sapphire filter to decrease the epithermal neutron component of the beam. To my understanding, NCrystal works for the thermal neutron range, and there is no tool that is capable of simulating what I need. My question is therefore, what exaclty would be necessary (cross-section data etc.) to use NCrystal for this purpose in Geant4. Is it a matter of providing NCrystal with the relevant data and using the existing code to load and apply it in the simulation?

I will be very grateful for your help!

Cheers,
Klaudiusz

[tkittel]

Hi Klaudiusz,

So if you are using NCrystal in Geant4 (or OpenMC for that matter), you should most likely already have what you want! When using NCrystal in Geant4, NCrystal only modifies the usual neutron scattering physics (specifically the hadElastic G4 process) for neutron energies below 5eV (currently this is a hardcoded threshold value, I can tell you where you can change it if you wish). So neutron scattering at E>5eV, or indeed neutron absorption physics at any energy, is still handled via the usual Geant4 processes. So I believe you should be good to go, in particular since neither Al or O seems to have resonances below 5eV. Assuming you plan to model your sapphire filter with Bragg diffraction turned off (as opposed to actually doing the more expensive single crystal modelling), here is how your sapphire filter cross sections look in Geant4 (v10.4 since I am stuck on an older Geant4 version on my machine for other reasons):

Btw., note that the Geant4 hooks are up for an overhaul (most likely expect something this fall), and currently have issues with running in MT mode or with parallel geometries.

[9klaud]

Hi Thomas,

Thank you very much for your answer and the details! It sounds great, I will run some simulations and verify the output against the literature.

In fact, I need the model to be as accurate as possible. Could you please suggest how I can turn on the Bragg scattering?

I am actually running all simulations in MT mode (i.e. due to a complex geometry that takes around 50 GB of RAM to load). However, initialisation and construction of the geometry (single world) are done on a single thread. Phase-space file read-out and the run itself is spread across several cores. Is it safe to use NCrystal in this scenario? Please let me know if you need more details.

[tkittel]

No, unfortunately I do not think the current Geant4-NCrystal bindings are safe in such a scenario.

For how to enable the Bragg scattering, you need to use the dir1+dir2+mos variables cfg variables (see here: https://github.com/mctools/ncrystal/wiki/CfgRefDoc), but do you actually know the orientation of the single crystal? It is normal to model single crystal filters with bragg=0, to emulate the situation where the incoming neutrons are not aligned with any plane.

Cheers,
Thomas

[9klaud]

I see, thank you for the confirmation.

I am waiting for information on which orientation is used at the physical beamline, but in the literature, 100 is the most common. Alternatively, the c-axis 001 can be used but is more expensive to produce. Nonetheless, I would test both orientations.

I have tested two things so far:

1. I ran a couple of quick simulations with the sapphire filter in free air on a Geant4 MT installation using both, G4MTRunManager with multiple cores and the single thread G4RunManager. Interestingly, I did not see a difference between the results. When I cout the sapphire crystal material definition in DetectorConstruction, I can see all the parameters pulled from the NCrystal data, so I assume the geometry is correctly reconstructed.

2. According to the literature, the calculated neutron transmission for a 7.5cm thick sapphire crystal (100) should be Tn < 8% for En > 1 eV and Tn > 85% for En < 0.02 eV ("the fluctuations in the transmission curve which are due to Bragg reflections are less than 5%"). However, the output from my simulation is Tn ~ 38% for En > 1 eV and Tn ~ 31% for En < 0.02 eV. This is when I look at the initial kinetic energy of the neutrons that reached a sensitive volume.
   When I look at the final neutron KE when entering a scoring volume, it looks like the crystal acts as a moderator rather than a filter and shifts the spectrum. Do you think that the usual Geant4 physics might be missing some interactions, i.e. multiphonon scattering for En > 1 eV? In any case, modelling Bragg reflections will probably not help in achieving a Tn of >85% for thermal neutrons, so as you suggested, I might pass on this one. I am not sure how I can improve the agreement though.

I was also wondering, is it possible to override the temperature of the crystal?

[tkittel]

Hi again.

So first of all, it seems to me from the numbers you quote above, that although you have created NCrystal-G4Materials in your Geant4 setup, you have not actually enabled the NCrystal process in the physics list. Have a look at the last lines in the example at https://github.com/mctools/ncrystal/blob/master/examples/ncrystal_example_g4sim.cc for how you must do this. Specifically, you must invoke G4NCrystal::installOnDemand() after runmanager initialize, and before the beamon. And once again, this will currently not work with Geant4 in MT mode.

Concerning the actual NCrystal settings for the sapphire filter, I thought it sounded interesting enough that I went and created a little example jupyter-lab notebook for you here (for general introduction to the notebooks, see the README here).

The point is that you can take advantage of the NCrystal API to investigate your physics and material configuration without having to actually launch something complicated like Geant4. Then, having verified which settings you want to use, you can then use those in an actual Geant4 job. It also helps you to pinpoint whether the issue is in NCrystal itself, or in the way it is used in Geant4.

Note that the notebook should be downloaded and run interactively, since the github preview does not really show anything interesting when it comes to an interactive widget. Here is a screen-shot as a teaser :-)

[9klaud]

Hello!

I apologise for not refreshing the topic for a while; it has been a hectic couple of weeks.

Thank you very much for creating the jupyter notebook! It's very handy indeed! I investigated the code and you were right - the NCrystal physics was not turned on despite the correct definition of the crystal structure.

I ran some more tests, and now the results seem to agree with the jupyter notebook example and the reference literature for neutron energies below 1 eV, which is great! However, for En>1 eV, neutron transmission remains at ~30% +/- 10%, instead of <8% (as per jupyter notebook and literature), but as you noted earlier, this comes down to Geant4 physics.

[tkittel]

Hi again!

Concerning the transmission rate for En>1eV, +-10% is a bit high, so perhaps it would make sense to run again with higher statistics? Also, note that when used in Geant4, NCrystal provides neutron scattering physics up to 5eV (this is currently a hardcoded threshold, independent of material), so for a neutron in the range 1eV..5eV, it should see the same scattering physics as in standalone NCrystal (e.g. the notebook). Of course, the absorption physics is 100% from Geant4, but I would be a bit surprised if that made a difference. If you wish to investigate the discrepancies further, it might make sense to use a fixed neutron energy (e.g. 4eV) for the comparisons. And more statistics.

Another thing worth noting is of course that modelling a sapphire filter as a high quality (sccutoff=0.2) single crystal might consume a lot of CPU power. One idea I have to improve this is to add a new special "single crystal filter" mode. So one would provide just the primary direction (e.g. 100 or 001), mosaicity, and energy resolution (e.g. delta-lambda/lambda) and NCrystal would calculate and cache the cross section curve during initialisation - and simply use that curve (+an assumption of isotropic scatterings) during simulations. Mind you, such a mode does not exist at the moment, but if there was enough interest I might be persuaded to look into it (although most likely not for ~6 months).

Cheers,
Thomas

[9klaud]

Hi Thomas,

In regards to the variation in the transmission rate of +/- 10%, I meant the transmission per energy bin - I used a neutron spectrum that extends up to a few MeV. I apologise for the confusion! It is far from the results in the literature though (<8%).

Nonetheless, I tested multiple monoenergetic beams from 1E-5 eV to 10 MeV. What I noticed is that the transmission rate follows the jupyter notebook widget and literature up to around 0.1 eV. From there, it looks like other Geant4 processes are taking over. Transmission rates were around:
76% for 0.1 eV;
52% for 0.3 eV (should already be ~35%);
46% for 0.5 eV (should be ~25%);
41% for 1 eV;
39% for 100eV;
40% for 1 keV;
41% for 10 keV;
26% for 1 MeV;
56% for 10 MeV.

I tested other parameters, i.e. sccutoff=0.2; dir1=@crys_hkl:1,0,0@lab:0,0,1; dir2=@crys_hkl:0,0,1@lab:1,1,0; mos=1arcmin. I did not notice any difference. It is interesting because the transmission for 1eV < En < 5 eV should yet drop. I also turned off the Thermal Scattering Libraries in Geant4, but it did not help either. As for the Geant4 physics lists, I use: G4DecayPhysics, G4RadioactiveDecayPhysics, G4HadronPhysicsQGSP_BIC_HP, G4HadronElasticPhysicsHP, and G4StoppingPhysics; 10 um cuts for neutrons.

Re modelling high quality crystal:
This is a great idea. I would expect the results to be closer to the literature at the moment. It is great below 0.1 eV, so fixing the agreement above this energy would be good enough for me to proceed further. My simulations already require a lot of computing power, so I am trying to keep it as simple as possible wherever I can. :)

I will try to dig deeper into this. Do you maybe have some suggestions on where to look in particular?

[tkittel]

So this is a bit puzzling. Off the top of my head, a few ideas/comments:

• While debugging this, use the "...;bragg=0" configuration, since it is maaaany orders of magnitudes faster CPU-wise than a true single crystal one.
• Since the actual cross sections are extremely energy dependent (cf. what happens if you disable the "crude_energy_smearing" option in the notebook), it might make sense to not use a completely monochromatic neutron gun - but perhaps smear the neutron energy with 1%.
• Did you verify that you have the same thickness (7.5cm) as in my notebook?
• The scattering physics should be completely provided by NCrystal up to 5eV, but in principle, absorption physics could depend on your physics list and geant4 list.
  • Try to use a prebuild reference physics list instead (like "QGSP_BIC_HP_EMZ" which I used above).
  • The difference could come down to the G4 version. Unfortunately, it is not trivial for me to currently test with a later version (it is the first thing I will fix after summer).
• Are you 100% sure all your neutrons are going through the filter, and that you don't have some of them taking another path to your scoring volume (i.e. try to reproduce with a very simple empty world with a single G4Box (slab) of filter material.

[9klaud]

Hi Thomas,

Thank you for the suggestions! I ran some tests, and I realised the scoring plane was bigger than the beam size, so it also registered higher energy neutrons scattered away from the beam path. Once I set an equal size, the results follow the notebook and literature.

It is a fantastic outcome, thank you so much for all your prompt responses, effort and help! I genuinely appreciate that. I will make sure to cite NCrystal in the paper.

Best Regards.
Klaudiusz