Hi all,
I am using a histogram energy spectrum input to gps in a macro file and I am scoring the energy spectrum 0.1cm from the source to check that the spectrum output matches the input. However, I find that G4 matches the mean energy well but underestimates the lower energy regions by up-to maybe 5%. I am trying to do detailed dosimetry calculations and it is vital I reproduce the energy spectrum exactly. Does anyone have any experience with this? This commands in the macro file are shown below:
/gps/ene/type Arb
/gps/hist/file hist_input.txt
/gps/hist/inter Exp
The input energy spectrum I am using comes from a EGSnrc simulation for a particular real life accelerator and has been thoroughly documented so I am sure that there are no issues with the energy spectrum. I also find that G4 essentially shifts the whole energy spectrum up by half a bin width, so I have to reduce the energy value in the hist file by half a bin width so that it matches the input spectrum from EGS. Here is the plot of the input and output spectra. EGSnrcSpectrum_G4RecreatedSpectrum_Comparison.pdf (14.9 KB)
Thanks for reading and for any help!
Michael
This is possibly a bug in GPS I found some years ago. My workaround was to convert the spectrum to an accumulated spectrum and use /gps/ene/diffspec false.
If this fixes it for you, please let us know. I will then submit a bug report.
What do you mean by an accumulated spectrum? I’ll try the command you mentioned and let you know!
Thanks!
Ah! Sometimes called “cumulative”, also “integral”.
I think GPS assumes it’s integrated forwards, so to speak. Say you have a normal “differential” distribution 1 2 4 6 5 3 2 1 0 0 at 10 energy points, the integral spectrum would be 24 23 19 13 8 5 3 1 0 0, i.e., at each energy E it’s the integral E>infinity.
Hi there, sorry, I’m not sure I follow still. I have tried summing your example at each point and am not getting the result you give, am I doing something wrong or is this a mistake?
Thanks!
Maybe I got my sums wrong - I did it quickly, hoped you would get the idea. What did you get? I also oversimplified - the simple sum of all intensities only works if the interval between the points is the same for all points. But (I try to be careful) the integral spectrum at energy E is the integral of the spectrum from E to infinity.
Attach your hist_input.txt. What are the units?
Hi, here is the spectrum file where the first column is the energy in MeV and the second is the fluence.fluenceSpectrumEGS.txt (3.6 KB)
Thanks for your help!
Mmm. Well it looks like a differential spectrum - the values are not always monotonically decreasing (as they would be for an integral spectrum). If it’s differential the units for fluence should be number per cm^2 per MeV, or some such. Can you confirm that? It’s not enough simply to say “fluence”. (If it’s integral it wouldn’t have the “per MeV”.)
If it’s differential, the values may not be the estimate of the fluence at a given energy. It may be a histogram - the value of the average over the next interval (or bin). You could try /gps/ene/type User.
What is the source of your data? You may be able to get hold of the equivalent integral fluence.
To calculate the integral fluence from a point sampled fluence at energy E(i) (say) do a step-wise integration, summing the mean of each interval, Sigma 0.5*(E(j)+E(j+1))*interval(j), j = 1,…,N. In your case, interval is always 0.1 MeV.
To calculate the integral fluence from a histogram, it would be Sigma E(j+1)*interval(j).
Sorry if this is all confusing. I recall being confused myself, and I’m not sure I’m totally making sense now.
Here are a couple of further comments:
Looking at my data (from long ago) I see I started at 0.1 MeV. Particles with energy below 0.1 MeV did not reach my detector. I wonder if GPS has trouble with the concept of a flux at 0. Try just leaving out the first data point.
All in all, integral spectra are more amenable to a precise definition and are used commonly in space physics.
Hi again, thanks for your reply! From what I understand it doesn’t really matter what the units of the spectrum are to GPS, as it is automatically being adjusted to whatever the required GPS input is, as long as the trend of the values is the same and describes the same spectrum shape. I’ll try to confirm the exact units of the spectrum I sent however. The source of my data is a phase space file from an EGSnrc simulation.
Yes, I am getting a bit confused by you’re equations above. I have still not managed to get the correct spectrum but I will keep trying with your suggestions and also try to get the integral spectrum from the phase space file and hopefully something will work. If I find the solution I will reply here with what I have done.
Thanks again for your help!
Hi again,
So I have tried what you have suggested with an input integral spectrum which looks like: EGSFluenceSpectrum_12MeV.pdf (12.7 KB)
I obtained the spectrum by summing the bin values from the bin E(i) to infinity.
Then using the commands:
/gps/ene/type Arb
/gps/hist/type energy
/gps/hist/file <integralSpectrumAbove>.txt
/gps/ene/diffspec false
/gps/hist/inter Exp
I obtain again the same underestimation of the the mid to low energy region of the spectrum: EGSnrcSpectrum_G4RecreatedSpectrum_Comparison.pdf (15.3 KB)
I also tried commenting out the 0 MeV entry or setting it to zero also and I get the same spectrum again.
Is this the type of spectrum you were referring to, or am I still missing something?
Thanks,
Michael
Hi Michael
Thanks for persisting with this. We need to get to the bottom of it.
Looking again at your fluence spectrum fluenceSpectrumEGS.txt, it looks like the abscissa values are the low-edge of the bin of the histogram. So integrating with x.cc (700 Bytes) gives integralFluenceSpectrumEGS.txt (3.6 KB). Does it match yours?
John
c++ x.cc
cat fluenceSpectrumEGS.txt | ./a.out > integralFluenceSpectrumEGS.txt
Hi John,
This spectrum you sent matches mine almost exactly. I am using the low-edge bin value because if I do not, the whole spectrum is shifted up in energy and does not match the input spectrum at all!
Cheers,
Michael
So the exact EGS output file is:EGSFluenceSpectrum_12MeV.txt (5.6 KB)
Where the first column is the mid point of the bin with units of MeV and the second column is the fluence with units fluence/MeV/incident particle.
Like I said before I have to reduce the energy value by half a bin width i.e. 0.05MeV in order for the spectrum to match, otherwise G4 somehow increases the energy of the spectrum.
Maybe I am getting something wrong at this stage by simply decreasing the value of the energy bin.
Michael
Hi again
This issue is troubling me, so I looked back at your original posting. Are we talking about electrons or photons?
You say you are “scoring the energy spectrum 0.1cm from the source”. A lot can happen in 0.1cm - an electron can lose energy, a photon can scatter. Better to pick up the energy as generated by GPS. Anywhere you have access to the event, for example:
void YourEventAction::BeginOfEventAction(const G4Event* event) {
G4double ke = event->GetPrimaryVertex()->GetPrimary()->GetKineticEnergy();
G4cout << "Kinetic energy: " << ke << G4endl;
GPS says, “/gps/hist/point Ehi Weight Specify one entry (with contents Weight ) in a histogram (where Ehi is the bin upper edge).” Your fluenceSpectrumEGS.txt looks like the energies are bin lower edge. You say you shifted the spectrum by half a bin width. Maybe that should be a full bin width.
Why do you choose Exp for /gps/hist/inter Exp?
That’s it for now.
John
Hi,
We are talking about electrons. I have also tried scoring the spectrum in vacuum as well as air and in both scenarios the spectra are almost exactly the same.
Yes, this is the most confusing part that the guide says to use upper bin edge, however when I add half a bin width (to give me the upper edge value of the full bin) the spectrum is shifted up in energy. If I reduce by a full bin width I will get negative energy value in the first bin.
I will try your suggestion inn the EventAction and see if I can obtain a more direct spectrum that way.
I have also tried all interpolation types such as Spline, Lin and Log and again I get back an almost identical spectrum each time.
Thanks and sorry that this is being such a nuisance,
Michael
Hi John,
Success! Using the high-edge bin energy value and cumulative spectrum, scoring kinE in EventAction and plotting the histogram gives a very good result. In python I can retrieve each bin value and in order to compare the shape with that of the input, I have simply done binValue/max(binValue) so that the maximum value will be 1 in the spectrum as shown:EGSnrcSpectrum_G4RecreatedSpectrum_Comparison.pdf (18.8 KB)
where the orange histogram is the density plot from EventAction output.
As you can see the input and output plots match almost exactly which is great news. This is clearly an issue with the accuracy of the flux or current scorers when used in such a way.
Is there any more details you would like about the result?
Thanks for all your time,
Michael
Hi Michael
That’s great. Phew!
Could you please do one more thing. Could you try going back to the differential spectrum - with a bin edge adjustment - and see if that gives satisfactory results? As I said in earlier exchanges I did at one time find some problems with that.
John
Hi John,
here is the final plot with both cumulative and differential spectra compared. EGSnrcSpectrum_G4RecreatedSpectrum_Comparison.pdf (21.7 KB)
It may be hard to see, but the differential one shows slightly lower values in the mid energy range, and the normalised spectra are close but not exactly identical.
Cheers,
Michael
OK. Great. Glad this is resolved all round.
I’m not surprised there are minor - very minor - differences, perhaps slightly different ways of dealing with the data. My judgement would be that these differences are absolutely negligible, since (a) there are bound to be modelling errors, and (b) most events, in your case, come from the 11.5 MeV peak.
Good trucking
John
Hi Michael,
I am interested in the energy fluence you mentioned, but I don’t know how to calculate the energy fluence, could you tell me?
Lee