How to keep track of neutral pion in decay chain?


I am developing Geant4 based simulation, and I would like to simulate proton beams impinging on graphite target, and to check how many neutral pions produced and thier energy and direction distribution.
But I cannot find neutral pion in the pdgID, and can find it in parentID but with very small step size and the corresponding pdgID is positron which is strange.

It is appreciated if someone tell me how to keep track of pi0s in Geant4. Reconstructing it from 2 gammas is also sufficient if possible and .

The PDG code for pi0 is 111. You should be able to find pi0’s as secondaries produced by the various hadronic processes. Since the pi0 lifetime is extremely short (~1e-16 s), it’s generally not tracked in Geant4; it’ll be created and take a single step to decay.

If you want to reconstruct it from the two decay gammas, you could either do that in your tracking action – look for two gammas from the same parent. Or you do it at analysis level, looping over the final-state gammas.

Thanks for the advice!
I ran 100 events and looked for gammas and check their parent, but pi0 does not exist … (i.e. pdgID == 22 && parentID == 111 does not exist)…
Do you have any idea about this?

I checked as
TTNNHit* aHit = static_cast<TTNNHit*>(i_TTNNHitsCollection->GetHit(i));
G4cout << "ParticleName: " << theEventID << " " << aHit->TrackID << " " << aHit->ParticleName << " " << aHit->pdgID << " " << aHit->parentID << " " << aHit->StepStatus << " " << fGeomBoundary << G4endl;

What physics list are you using? Do you have appropriate hadronic physics? What’s your beam energy?

It’s possible that you’re either below threshold for pi0 production, or you just didn’t run enough events.

I tried FTFP_BERT and CHIPS, but situation looks similar.
How can I check included hadronic physics?
Proton beam energy is 30 GeV so it is enough to produce pi0.
hmm, maybe I will try more event as 10k.

Both of those are hadronic-based physics lists (FTF is for very high energies, BERT (“Bertini”) is a nuclear cascade model suitable for about 20 GeV or below. You’re quite right that a 30 GeV proton beam on carbon (CMS energy is about 19 GeV) is plenty high enough for pi0 (and K0 and rho and . . .) production.

I don’t have an intuition for the production rate, though.


I increased the number of events from 100 to 1000, then pi0s appear in pdgID.
(But still no pdgID == 22 && parentID == 111, maybe the step size 1 mm is too long and next process happens during the 1 mm.)
Another strange thing is that there are multiple pi0s in an event, and they have same coordinates, but different kinetic energy, momentum and track ID.
Do you have any idea about this?

The coordinates would be the same if they all came from the same step (post-step point). Geant4 doesn’t track positions within the nucleus of an atom. Do you store the track number and step number of the parent when you collect the secondaries?

Most hadronic interactions at the 10-20 GeV level will produce lots of secondaries, so I’m not surprised you see multiple pi0s. There are probably a bunch of pi+ and pi- also.

Finally, I didn’t read your if-condition clearly. It’s probably wrong. “parentID” is almost certainly the track number of the parent, not the parent’s PDG code, unless you wrote your SD to collect the latter specially.

Thank you for the explanation (especially parentID that I completely misunderstood…)
And yes, now I come to think multiple pi0s maybe produced from a single identical process. (though at maximum 10 pi0s are produced at the same position…)
So parentID is track number of the parent, right?
And I did not keep parent step number information, so I will try to include it.

Yes, parentID is the track number of the parent (you can get it from G4Track::GetParentID()). Note that once the parent track is finished, it is deleted so you can’t use that ID to go get further information. If you want to store the parent PDG code or anything else, you’ll have to do it with a buffer in your SD or SteppingAction.

Here a macro for example TestEm1 and its printout with /tracking/verbose 2
An eta of 1 eV decays in 3 pi0
tomoya.mac.txt (525 Bytes)
tomoya.out.txt (7.6 KB)


Thanks! First I was surprised that 10 pi0s are generated at the same coordinate, but maybe it was just an extreme case.
Now I think my simulation works correctly (i.e. without any unexpected and bug).
Thank you for your help!