Hi all,
I’m using Geant4-v11 and it seems to me that gammas coming out of positron annihilation are always perfectly collinear, is that correct? I’m a bit puzzled, because, when I look into the G4eplusAnnihilation class, it seems that the G4eeToTwoGammaModel is initialized, which takes into account a possible in-flight annihilation of the positron. However, that model class is never called when I run my code.
How can I simulate correctly the non-collinearity of gammas?
If I understand correctly from the output, in this case the simulated positron annihilates with zero momentum, therefore at rest, not in-flight, and the two produced gammas have not 511-keV energy, as I would expect.
I’m a bit puzzled by this: if the positron annihilates at rest, the two gammas should be emitted back-to-back and with 511 keV each, right?
What I would like to simulate is positrons annihilating in-flight, when it is the case, therefore the gammas should be emitted slightly non-collinear (in that case).
No. It is a annihilation in flight.
At end of step 11, the positron as a kinetic energy of 21.15 MeV
Along the step 12, it loose 2.564 MeV, so that its kinetic energy at end of step is 21.15 - 2.564 = 18.586 MeV.
So that the total energy available for the annihilation is 18.586 + 2*0.511 = 19.608 MeV
This is the total energy of the 2 gamma created.
I see, I was misled by the fact that the kinetic energy of the positron in step 12 is listed as 0… Is that incorrect, then, right?
However, I’m using the G4EmStandardPhysics_option4 physics list, and the gammas are always generated perfectly collinear. I’ve checked that for positrons only the AtRestDoIt method of the G4eplusAnnihilation process is called (no PostStepDoIt method of the base class G4VEmProcess seems to be called), and in that process gammas are generated back to back. Which process is the responsible for generating the gammas in this example?
Should I use a different physics list, maybe? I’ve looked in the PhysListEmStandard class, which is used in this example, if I understand correctly, and I don’t see anything special for positrons.
Thanks again! I managed to run the TestEm1 example and I find indeed the same result as yours: some positrons annihilates at rest, and some of them annihilates in-flight. I simulated a block of aluminum in my application and I find the same behaviour. However, when I inspect the output file, I find the attached distribution: in most cases, the gammas are perfectly collinear, while in the rest they have very large angles (it is not a smooth distribution, but I don’t know what to expect for this simulation).
I turned back to my applications, which is PET scanner. It is measured that the angle between the two gammas for typical radioisotopes (such as F-18) in body tissues is a gaussian distribution centred in 0 with a FWHM spread of 0.5 degrees (see, for instance, DOI: 10.1007/978-3-540-36841-0_411). I repeated my simulation using water (the density is similar to a body tissue, therefore it is commonly used in PET simulations) and an energy for the positron of a few MeV, which is the Q-value of many of these radioisotopes. Instead of a gaussian with 0.5 degrees spread (or something similar) I still find a distribution where in most of the cases the gammas are emitted collinear, and sometimes (in 7% of the events) they are emitted at very large angles.
It is certainly not the behaviour I would expect from physics, but maybe I’m overlooking something.
Thanks everybody for your help!
I am not sure to understand what you say.
Did you have a look to examples/basic/B3/B3a ?
About F18, run example rdecay01 (in radioactivedecay subdirectory) with attached macro.
Yes, I mean opposite directions, sorry! The plot represents 180 - angle between gammas. I plotted that way because it was easier to check my expected distribution.
I ran the example you mentioned and I find the output attached. The positron seems to end up in the World volume and no gammas are produced. Looking at the geometry, it seems this is just a big volume of air, so maybe the positron exits the volume always. output_f18.txt (30.3 KB)
In example rdecay01, only radioactive decay is registered; no other interactions.
My purpose was to show you the energy spectrum of emitted e+ : max= 627 keV, mean= 250 keV
I see. To be more accurate, I have simulated directly F-18 in water and the angle between the two gammas is not what expected. I’m wondering if I need to write my own annihilation process to produce the correct momenta of the gammas or if there is something missing in my simulation.
Thanks a lot
I do not know what you got, nor what you expected …
I tried to repeat your exercice with TestEm1 : add an histogram in HistoManager, and the following code in SteppingAction :
//plot (e+ e-) annihilation
//
if ((process)&&(process->GetProcessName() == “annihil”))
{ const std::vector<const G4Track*>* secondaries
= aStep->GetSecondaryInCurrentStep();
if ((*secondaries).size() == 2) {
G4ThreeVector gamma1 = (*secondaries)[0]->GetMomentum();
G4ThreeVector gamma2 = (*secondaries)[1]->GetMomentum();
G4double alpha = gamma1.angle(gamma2);
analysisManager->FillH1(8,alpha);
}
}
I attach macro and plot. Indeed, this is the final state generated by the annihilation process. F18.mac.txt (457 Bytes)
What I simulated was F-18 decays in water.
What I got is a distribution similar to the one we both find in example TestEm1.
What I expected is a gaussian centred in 180 degrees with a FWHM of 0.5 degrees (or a bit less since the value for pure water is lower than for body tissues), according, for instance, to Limit of Spatial Resolution in FDG-PET due to Annihilation Photon Non-Collinearity | SpringerLink.
I’m sorry to insist, but I would like to ask Geant4 developers if this is a known issue and if I need to write my own process or there is something I’m missing here.
You can’t get a distribution centered at 180 degrees in theta. 0 deg and 180 deg are special values in angle. If you measure the opening angle between tracks, 180 deg is hard maximum.
You can’t get a distribution centered at 180 degrees in theta. 0 deg and 180 deg are special values in angle. If you measure the opening angle between tracks, 180 deg is hard maximum.
Yes, of course, I didn’t express myself correctly. What I meant was that the expected distribution is a half gaussian, with maximum in 180 degrees and sigma 0.21 degrees.
But that’s not what happens. You won’t see half a Gaussian with hard edge. What you should see is a vaguely “Maxwellian” shape, with a minimum at 180 deg, a fast rise to a peak at less than 180 deg, then a longer tail outside that.
You can try it yourself with a simple random generator. Generate a Gaussian distribution of signed values X centered on zero (mean=0, sigma=1). Then make a plot of abs(X) from 0 to 5. You’ll see what I’m talking about.
I’m not sure what you mean… I’m claiming that this is the expected distribution according to what I found in literature and in several courses of nuclear medicine. See, for instance, the conclusions of the article I mentioned above: “The non-collinearity of electron-positron annihilation radiation from in vivo 18FDG was evaluated for the first time by the conversion from the photopeak spectrum into the angular deviation. The distribution was well described by a Gaussian function centered at 0 ̊ with the FWHM of 0.54 ̊”.
What they do is to convert the spread in the energy to a spread in the angle, thanks to the relationship: theta = 2 Delta E / mc**2.
