I’m trying to characterize some neutron interactions in a water volume (in order to compare the effect of different physics list and of the neutron energy) and I’m having some difficulties to understand some Geant4 output. Maybe some experts of hadron simulation could answer some of my questions?
My detector geometry is just a very large volume of pure-water water (with H and O natural isotopic distributions). I got the expected H neutron capture interaction and some O neutron capture interaction:
i.e. O or H nuclei whose creator process is a “nCapture”
And I sometimes got what I believe to be absorption process: the neutron track ends with an “inelastic” process and no neutron in the secondaries particles. I noticed something confusing when looking at the secondary particles in this case:
I most of the cases there is one nuclei (13C, 15N, 16N, etc.) and the interaction is quite simple to understand.
However in some cases there are multiple nuclei in secondaries (very often a O and a C nuclei but not always), all of them having as parent track the neutron. Here it’s a bit more difficult to understand what’s going on. It’s not possible the neutron break a bigger nuclei to produce them, since there a no bigger nucleus than O in the simulation. The only reasonable explanation I could think about was that these “additional” nuclei in the secondaries are due to nuclear recoils occurring right before the neutron absorption. However I was expecting Geant4 to have one step per interaction, and so have one step for each nuclear recoil, it seems to not be the case, would you think this interpretation is correct? And is it the expected behavior? Can it be change with some parameters?
Giving this situation it’s sometimes quite hard to know with which type of nuclei the neutron did the main interaction (H? O? If O, which one?). Is there a way to know with which target nucleus did the neutron interact? I look in the different method of G4Track and G4Step but it seems there is no such function but I may have overlook.
“The only reasonable explanation I could think about was that these “additional” nuclei in the secondaries are due to nuclear recoils occurring right before the neutron absorption.”
This is not possible in Geant4.
You should tell in which G4 version and with which physics list or physics constructor you see these phenomena. If it is ParticleHP, which environment variables do you activate ?
Normaly, you would be able to reproduce them easily with example Hadr03 and to post the output.
Thank you for your answers. I am glad to know this is not the expected behavior.
The version of Geant4 I was using for this simulation is rather old (9.6p04), the physics list, for inelastic processes is the following:
G4NeutronInelasticProcess* theInelasticProcess =
new G4NeutronInelasticProcess("inelastic");
G4NeutronHPorLEInelastic* theInelasticModel = new G4NeutronHPorLEInelastic();
theInelasticProcess->RegisterMe(theInelasticModel);
G4NeutronHPInelasticData * theNeutronData1 =
new G4NeutronHPInelasticData;
theInelasticProcess->AddDataSet(theNeutronData1);
G4CascadeInterface* bertiniModel = new G4CascadeInterface();
bertiniModel->SetMinEnergy(20*MeV);
theInelasticProcess->RegisterMe(bertiniModel);
The following variable is activated:
setenv G4NEUTRONHP_USE_ONLY_PHOTONEVAPORATION 1
I tested Hadr03 example with geant4.10.03.p03 or with geant4.9.6.p04 using HadronPhysicsQGSP_BERT_HP or G4HadronPhysicsQGSP_BIC_HP and I couldn’t find the behavior I was observing in my original code. Likely the issue comes from my physics list or from some setting regarding how the steps are managed (but I don’t remember changing anything regarding this).
Anyway, looking in Hadr03 example I find the answer to my second question, we can know what was the target of the interaction by using something like this:
G4VProcess* process = const_cast<G4VProcess*>(theStep->GetPostStepPoint()->GetProcessDefinedStep());
G4HadronicProcess* hproc = dynamic_cast<G4HadronicProcess*>(process);
const G4Isotope* target = NULL;
if (hproc) target = hproc->GetTargetIsotope();
G4String targetName = "XXXX";
if (target) targetName = target->GetName();