I have a simulation where I am shooting neutrons at a block of quartz (SiO2). I expect a neutron spallation reaction to happen such that n + 28Si → 26Al + 2n + p + e- (production of 26Al which is used for cosmogenic nuclide dating). I need to know how many 26Al atoms were produced in my simulation. Is there any way I can call a function to just get the number of certain atoms in a physical volume? I looked for such a function but could not find anything.
The way I am attempting to do it now is, just detect the 2n + p + e- from reaction in a detector. The issue with that is, I am pretty sure I am also detecting other reactions which also produce 2n + p + e- which is throwing off my 26Al number.
If you set the “proton” production cut to zero (or very small, like 1e-6*mm), that will apply to all hadrons, including nuclei. In that case, you will get the 26Al (or whatever else) nucleus created as a secondary track. You can then count it directly in your SD, or SteppingAction, or StackingAction, whatever you prefer.
In a solid material, the secondary nucleus will typically take just one step before stopping, so this won’t affect your CPU time significantly.
I am new to Geant so it took me a bit to figure out why your solution would work and how to implement it. I just want to check if you think this solution should work. This is my UserSteppingAction class I made to implement your suggestion.
I just have a counter up in my run action class to keep track of the number of 26Al I saw thought the entire run.
The reason I ask if you think this would work is, I am getting results way higher than I expect (I am counting much more 26Al being produced from 28Si than I expect. Getting 532 even though I expect 6-7). I suspected with this code I might for some reason be multiple counting the same atom. Another possibility I thought about, if a neutron spallates 28Si to 26Al, the 26Al has a decent chance of first ending up in its non ground state. With this code, would I also be counting each de-excitation of the 26Al all the way to its ground state. So for example, lets say when the 26Al was produced it was in a state around 2000KeV above ground. The 26Al nucleus then de-excites to a state around 400KeV above ground. It then de-excites again to the ground state. With this code, would I have just counted that 26Al 3 times (one for the initial production, one for the first de-excitation, one for the second de-excitation).
Is this all true? If so it makes my problem easy(printing out the tracks; I am detecting the same track IDs many times over). I am a bit skeptical the tracks are truly unique across the whole run of the simulation regardless if the original track was killed. Printing out my track numbers I am getting numbers in the singles to 10s range. In my simulation I am sending 10 million neutrons, so if each of those neutrons and their secondaries each have a unique track ID, I would expect to get track IDs in the millions range. I will say, making it such that I cannot detect a track ID more than once, my number of 26Al production is now 19 which is much closer to what I expect.
The track IDs are not unique across all events; they are guaranteed to be unique (and increment starting with 1) within each event. So if you save track information to your N-tuple, I strongly recommend saving both the event ID and the track ID.
What you did in your SteppingAction is exactly right, as far as it goes. The first-step test guarantees you only single-count the track within an event. However, your point about excitation is right on. The 26Al will almost certainly have been produced in an excited state, which then has to be passed (via tracking) to the de-excitation process. G4 won’t treat the intermediate states as separate tracks (all the de-excitation to ground happens in one step, with multiple gamma emission). So that might get you two counts per track.
One way to test that would be to add a debugging message – in your if block, print out the particle (particle->DumpTable(), I think) and see what you get.
You have been a huge help. I have one more question regarding GEANT4 but it requires I very briefly explain what I am trying for my project I am working on.
I am trying to simulate the production of 26Al from 28Si in quartz due to neutron spallation. The reason I am trying to do this is, 26Al is used in geological dating and its production rate needs to be well known. Therefore, I am trying to do that via simulation (95% of 26Al is produced from neutron spallation). The way I am testing if this works is, I am taking a real place on earth where 26Al production rates are well known, simulating the neutron flux and energy density at that location incident on quartz (I just have a 1cm^3 quartz block in GEANT) and seeing how many 26Al are produced.
I am simulating 1 year’s worth of neutrons in GEANT4. From real measurements, I should see 209 26Al atoms produced in 1 year’s worth of neutrons, however, in my simulation I see about 43,000. The most obvious thing that could be happening is I am somehow multi counting my 26Al atoms. However, when I go and print out the event and track ID of my 26Al detection, I see the eventIDs differ by like hundreds or thousands so I thought it unlikely I was multi counting (else I would see my detection happening in the same eventID). Another thing I thought was odd was, I was running DumpTable for each 26Al and almost all of them were in their ground state.
The fact almost all of them are in their ground state makes me scared I am somehow multi-counting because I would expect many more of them to be in excited states (which I did confirm a very small amount of them were).
My question is, could I somehow be detecting the same 26Al between different events? I was assuming the way Geant4 works with events is
1: sends in particle from the particle gun
2: all interactions happen and secondary particles created
3: all interactions will eventually happen until all the particle tracks are dead
4: increment the eventID and send the next particle from the particle gun and repeat
Is this accurate or is it likely I am somehow detecting the same 26Al between the events?
Just some additional detail that might be critical. My physics list is QGSP_BERT. Also I never did set the “proton” production cut to zero because this method just seemed to start working without that.
You’ve described your simulation well. You are absolutely not detecting “the same 26Al between events.” Geant4 deletes all the tracks and starts over with every event. That’s extremely well validated.
I suspect that the problem is your normalization. You write, “I am simulating 1 year’s worth of neutrons in GEANT4.” How many events is that? And are you generating the neutrons as a beam, such that every one of them is incident on your silicon target? In nature, the neutrons are typically isotropic, so only a fraction of them would hit your target.
Okay, as I expected it is most likely somewhere the physics of the situation is wrong (My understanding of what I should be doing is flawed somewhere). Since this is a physics question I will close the topic since I have professors at my university I can further talk with. This project is just a senior project and I have already hit my goal of where I was supposed to get to so it is okay to leave it in this “unfinished” state. Thank you so much for your help.
But to answer your question, the way I am simulating the neutrons is with a GPS. I found a paper which gave me the neutron flux as a function of energy. That equation I use for my neutron energy distribution. The way I determine one year’s worth of neutrons is that I integrate the flux as a function of energy equation over the energy range I plan to simulate and get a total neutron flux out. I can then use that total flux to see how many neutrons occur over one year (the number comes to around 81 million neutrons). All my neutrons are incident directly on my target and all hit the target. However, I figured this is fine because regardless if the neutrons come in isotropic, the target should still see the same neutron flux (which I have gotten from a paper). One thing that could be wrong is, out in nature, the 26Al is produced inside small quartz crystals inside rock (in particular I am looking at the SP lava flow in Arizona). My sample in Geant4 is just pure quartz. My pure quartz sample doesn’t take into account the effect of all the other rock and stuff that do exist in nature.
That’s an excellent description! Yes, I’ll bet if you surrounded your quartz target with a volume of rock, you’d find the 26Al production rate would go down. All those other nuclei can scatter and capture neutrons as well, reducing the “effective flux” to your target. We do similar studies for backgrounds and activation in my experiment (a dark matter search down in a mine), and getting all the intervening materials right, along with the flux, is the hardest part.
