Good afternoon,
if my understanding is correct, for a low energy neutron the choice of the isotope and the reaction channel is done in
source/processes/hadronic/models/particle_hp/src/G4ParticleHPChannelList.cc: 165-176
How can I retrieve the result of this choice during transport? The only information that I have is hInelastic, type 4 and subtype 121. Can I retrieve the MT channel? Thanks and best regards,
The TMX Coder
Dear TMX-Coder,
One way to try to access to the isotope where the reaction took place is to have a look to the “G4ParticleHPReactionWhiteBoard” class. It takes the Z, A and M (isomeric state) of the nucleus in G4ParticleHPChannel.cc, lines ~230-235.
I don’t know how to get the MT channel, but probably it is not easy.
Regards,
Emilio
Hello Emilio,
thanks a lot for your answer. Here is my understanding of the situation. Effectively around line 243 and successive the target isotope is chosen. We enter this code only when the step is limited by a neutron interaction. If it is fission, elastic or capture, the matter is handled differently, if it is inelastic, this is done (AFAIU) in the following way:
in G4ParticleHPChannelList
there is a first double for where the data for each isotope and all the inelastic channels are summed. Then the isotope is chosen from the weighted (i.e. macroscopic) total inelastic x-sec, obtained via the sum of individual inelastic cross sections via G4ParticleHPChannel. It is not entirely clear why this is necessary, since each isotope has its own inelastic x-section already stored in the G4ParticleHPInelasticData, which could be used to select the isotope. Even better, one could select the isotope at the time of checking the path-length and remember the selection.
Then, once the isotope is chosen, the channel is selected with the same procedure. Note that storing all the channels x-sec for all the isotopes when selecting the isotope would avoid the costly procedure to recalculate the x-section at the price of few tens of floats. There are fast-paths to not recalculate the x-sections if all stays the same, but still a lot of overhead would be avoided. The fact that the first time the x-sec are macroscopic would not alter the selection process if properly normalized, since there will be the same constant for all channels of one isotope.
The way the x-sec os calculated for each isotope (in principle)
Sample the (biased) thermal velocity of the nucleus according to the Maxwell distribution;
Boost the neutron in the reference frame of the nucleus
Calculate the x-section at 0 degrees;
Sum all x-sections;
I see at least three problems with this code.
Minor. Sampling of the Maxwell distribution is done sampling px, py and pz and using rejection. It would have probably been better to sample the velocity and then rotate randomly.
Major. For some reason that I do not (yet) understand, the choice of the nucleus is “biased” toward nuclei moving in the same direction of the neutron, reducing the relative energy. I see no reason for this. The choice should be biased indeed, but by the cross section under consideration. If the energy of the neutron is just “below” a resonance, the neutron will “prefer” nuclei moving toward it, while if the energy of the neutron is just “above” a resonance, the neutron will indeed prefer nuclei moving “along” or away with it. There is no good way to do this in fact. Unless I am missing something…
The problem is that the information on the chosen channel is lost in G4ParticleHPChannelList. One way to make it “go back” would be to store it in G4HadFinalState, which is returned, but not really readily accessible. Another way would be to store it (in some way) into G4HadronicProcess, like the target nucleus.
To access the channel by channel x-sec is perhaps less cumbersome. G4HPChannelList is a static class. If one would add a getter like
G4ParticleHPChannel ** HPChannel() const {return theChannels;}
it would be possible to access the x-secs. Let me know what you think. Really sorry for the long mail!
The Transmutex Coder
Transmutex SA
Hello Emilio,
I understand that my rant was a bit difficult to answer to. But now I have a very precise question. In lines 249-252 of G4ParticleHPChannel I have
xsec[i] = theIsotopeWiseData[i].GetXsec(aThermalE.GetThermalEnergy(theTrack,
theFinalStates[i]->GetN(),
theFinalStates[i]->GetZ(),
theTrack.GetMaterial()->GetTemperature()));
I think there is a problem since the N and Z should be the ones of the target isotope and not of the resulting final state of the reaction. If my understanding is right, this is a serious bug. If not, I would be pleased if you could let me know what am I missing. Best regards,
Dear TMX-coder,
I have not written this code and I don’t know all the details. There could be bugs in it, of course, but I don’t think that this is one of them: the impact in the results of a simulation would be quite big.
Concerning your previous post, my answer is more or less the same. I don’t know all the details of the code and it would take me a long time to understand why it was written one way or another.
I’m sorry of not being of more help,
Regards,
Emilio