Hi,
is there any feature to receive the statistical error resp. uncertainty of a simulation. Afaik Monte carlo simulation usually provide quantities like these(?), is there also something available for G4?
My actual purpose is to determine how many photons I should generate to get stable results.
The answer depends on how you are getting your results. If it is via scoring , Geant4 provides stat. error, if you make an histogram, for example in ROOT, ROOT provides error bars in each bin
Thank you, yes, I’m making histograms via ROOT. How do I get an error bar, there?
I’m reading out a tree and creating a tree quite directly.
If you fill a histogram from a branch, then ROOT will compute symmetric error bars for you. You can have them drawn by specifying one of the “E” options to Draw().
If you do your own analysis, then you get to do stuff like compute sum of weights, sum of squared weights, etc. Woo hoo…
Thanks a lot, that worked quite flawlessly.
For my understanding: What exactly do these error bars represent?
For an unweighted histogram, they are sqrt(N) (“RMS”, “one sigma”) for each bin. For weights, they are sqrt(sum-weights) for each bin. See the ROOT documentation for more details, and for options to display error-on-the-mean, or asymmetric (Poissonian) error bars.
As usually, thanks a lot for all your help and inputs! The question is cleared but I’m still unsure of how many photons (or particles, generally spoken) I’ve to generate. How can I determine this?
Is there a thumb of rule for G4 or purposes? Can I use some calculated values?
Geant4 generates events, not “activity,” so what you’d need to do is understand first what your source activity is, and what time interval you want to simulate. The product of those two is the number of events.
Having said that, the next question is what uncertainty do you want (1%, 0.1%, 10% whatever…)? The relative uncertainty goes like 1/sqrt(N), so you can use that to estimate the number of events you need. Given the number of events, you can rescale to time interval or activity as needed.
Thank you! I assume by source activity you mean how many photons the source is emitting? I’m controlling this parameter, right? 
The time interval refers to the flight time of each photon? At the moment the time interval is set such the photons are killed when entering a detector the world boundary. How do I put these into numbers?
The level uncertainty is a good question, I’d make this dependent on the effort to reach a certain one. But in general I guess something between 0.5 - 5 % would be good.
“Activity” describes something in real life, for example “a 5 mCi Am-Be source”; it is typically in units of number of decays or number of emissions per second. In the simulation, you control the total number of events generated; you might need or want to translate that into a meaningful real-life source quantity like activity or exposure time or something.
Ok, got that. But what is the difference between events and activities, here?
E.g., when running the simulation I set /run/beamOn 10000 to create 10,000 photons, or? So I’d assume this is my activity?
Just to make sure: I’m not using decaying material but photons right from the start.
this is the product of activity [events / s] and time interval [s]
think of a very low activity source, then your simulated time interval would be very long (until that number of decays would happen) or a very high activity source, then the time interval is short.
That is fine, because you know how many photons you get per decay.
this is the product of activity [events / s] and time interval [s]
Thank you! How do I find out what activity and time interval I’m using? My code that, from my pov, produces the photons is just:
void MyPrimaryGenerator::GeneratePrimaries(G4Event *anEvent){
G4ThreeVector pos1(-53.66*cm, -30*cm - 1.8*mm, 0*mm);
fParticleGun->SetParticlePosition(pos1);
G4double phi = 78*deg;
G4double a = phi*G4UniformRand() + 90*deg - phi;
G4double b = 2*deg*G4UniformRand() - 1*deg;
G4ThreeVector v(std::cos(a), std::sin(a), b);
fParticleGun->SetParticleMomentumDirection(v);
G4double energy = InverseCumul();
fParticleGun->SetParticleEnergy(energy);
fParticleGun->GeneratePrimaryVertex(anEvent);
}
while InverseCumul() is a certain function to take into account a provided spectrum.
Hello,
Geant4 has no time interval in which your 10,000 particles are created.
Geant4 simulates particle per particle but they all start at time 0.
You have to introduce a time component yourself then do the math.
Let’s say you create 10,000 photons.
You are free to say these events happen in 1 second → then you have an activity of 10,000 Bq.
You are also free to say these events happen in 10 seconds → then you have an activity of 1,000 Bq (assuming constant activity) and your results are what you would measure after 10 seconds.
Best regards
Max
Thank you! Atm I don’t need this time information but I wonder how to implement it. I would simply let a variable count up but this is probably up to the computing system how much time it needs. I guess it is somehow possible to connect to the computer clock?
However: When G4 is not using activity but events, what is it finally doing (or using)?
My program stops after 10,000 produced photons (?!).
Hello again,
Now this sounds like you want to know the time your simulation takes to execute, which has nothing to do with converting “produced particles” to meaningful physical units like “activity”.
Let’s say you shot a proton into something. There are many interactions it can make e.g. (p,n), (p,2n), (p,d) …
All these interactions have different probabilities (cross-sections). You could roll a dice to decide which interaction takes place and do the calculation by hand but a single particle can do 1000’s of interactions. G4 does the dice roll for you and these calculations so you are able to calculate the output of millions of protons in a short time.
Regards
Max
Ok, thank you!
One last question: When G4 is taking into account such events instead of …the number of generated particles(?): Why does my simulation stop after the generation of all set photons?
When you start your run with /run/beamOn 10000 G4 creates one particle and proceeds it and all secondary particles created by it until the end (either the come to rest or they leave the world volume) before starting the next particle.
Your simulation stops after 10000 photons because everything is done.
And these 10,000 photons are 10,000 events, then?
An “ event ” in GEANT4 includes the results of one primary particle (as well as all its secondary particles) travelling through the geometry.