GetTotalEnergyDeposit() from step gives me a zero when filtered for particle type and whether the particle is within a volume.
However when I GetKineticEnergy() in pre and post step point and subtract the values and feed them into event action where they are summed it looks like I am able to obtain the energy spectrum of the particle of interest in the volume.
Why is GetTotalEnergyDeposit always 0?
The energy deposited in any detector is only from charged particles. Neutrons and gammas “indirectly” deposit energy through collisions/interactions with charged particles like nuclei or electrons. If you summed the total deposited energy from a neutral particle for all of its secondaries it would give you the expected values.
So doing it via the delta of the kinetic energy I described is correct for neutrons? To me it seemed like GetTotalEnergyDeposit accounts for all processes.
For neutrons the primary means by which they transfer energy to charged particles is by nuclear collisions/reactions. The process goes like this:
If you filter on neutrons you will lose the ability to account for deposited energy from (4). It will be appear missing. With charged particles, most of the energy loss is actually in (1) which would be part of total energy deposited. Typically sensitive detectors record all energy deposited by all particles to avoid that problem. You could also set the range cuts for electrons, protons, positrons, and gammas to a kilometer and that would push nearly every secondary into (3) if you wanted but this approach is very inflexible.
I think I am happy to get the energy of the neutrons I am interested in from the GetKineticEnergy function, though what you write isn’t immediately obvious in the manual. I should take another look. Thanks for your help.
I have another question, that I am struggling with. When I obtain the deposited energy in the volume im interested in from each primary particle, I get the same number of entries in my histogram as primary events due to no primaries being fully absorbed prior to reaching the physical volume of interest.
However when I try to obtain the energy of those particles in tracking action only once and on entry to the Physical Volume I am interested in, I only get a number of entries equal to the number of reaction events. I am filtering on the primary particle type and volume pointer. Shouldn’t this give me the energy of the primary on entry? Instead it seems to give me the energy of the primaries that ended up causing a reaction, which is a bit useful but not quite what I want. Any help with solving this logic error would be most welcome.
Cheers
What experiment are you trying to simulate? It will help to answer your questions since the answer might depend on what you are doing.
Expected.
You cannot use TrackingAction for this. It is only called at particle creation and particle destruction including AtRest processes. The reactions happen in your volume so that is why they are likely passing your logic check.
You will have to do this in the SteppingAction. Some variation of:
void SteppingAction::UserSteppingAction(const G4Step* step){
G4LogicalVolume* volume = step->GetPostStepPoint()->GetTouchableHandle()->GetVolume()->GetLogicalVolume();
const G4int particleID = step->GetTrack()->GetTrackID();
// Your code here...
if ((step->GetPostStepPoint()->GetStepStatus()==fGeomBoundary) && (volume==fVolume)){
// fVolume is a pointer to the logical volume you care about, passed from detector construction
// First statement checks that you are about to enter a new volume, second checks that you are entering the volume you care about
}
}
You can get the energy there then use the PostUserTrackingAction to find the energy when the particle is killed either from reactions or leaving the world volume. You could also throw in an added check that its your primary. A far more better way to do that is to just check the TrackID which is zero for primaries and no other particles. I have included it above for the stepping action.
A secondary neutron field that is generated by a reaction and to which I want to expose a sample. I have split the geometry into a simple case to check I am able to record the data I want from the theoretically produced neutrons and a more realistic geometry that generates said neutrons and record the reaction product I am ultimately interested in. That is why I want the spectrum of the neutrons entering the sample and ultimately the intensity as a function of the position and the effective spectrum in the sample (probably as a function of sample depth).
This is quite a good insight. Thank you.
I am also trying to develop my ability in geant4 so I do want to get to grips with the ins and outs of the toolkit and its nicer to have someone to talk to (apologies for any lazy questions though I am trying to avoid it).
I will try this and update.
This is exactly the use case for a SensitiveDetector object. You write it yourself, so that you can record exactly what you want to record, with whatever if statements you need. Almost every example has one; you can find an example and pick it apart with a find pipe like this:
$ cd <my-geant4-area>/examples
$ find . -name '*.hh" | xargs grep -l 'public G4VSensitiveDetector'
That will show you every single SD class included in every example. Pick the one you like!
In your SD, you can record the energy on entry to your volume by testing the G4Step preStepPoint for fGeomBoundary status, and record the kinetic energy odf the preStepPoint. You can record the energy deposit during the step by calling G4Step::GetTotalEnergyDeposit() (which is generally zero for neutrons, but might not be if you have production cuts!). You can collect the energy transferred to secondaries by looping over the secondaries attached to the G4Step and recording their kinetic energies.
Thank you for the tip 
. I also use grep -iRnw "string" ~/GEANT4/examples to look for relevant snippets along with telescope in vim. I should probably collect and categorise them at some point.
This is the task I will work on right after. I have found for now B4d and B5 along with RE01 that Makoto mentioned in his talk on scoring in the advanced course.
I came accross this in NeutronSource but I didn’t quite understand how to use G4StepStatus. I have simply been checking the current volume by GetStep()->GetPostStepPoint()->GetTouchableHandle()->GetVolume(); and comparing it with the pointer of the volume that a function in DetectorConstruction returns to the class I am using. I am not sure if this is perhaps uncouth but I will definitely try the status method as it seems abit more elegant!
I have set a fairly low 1 meV SetMinEnergy for the neutrons as the object is to produce a thermal field, but while in the sample I don’t expect more neutrons to be produced (some are but at high energy ). Either way this tripped me up already and I was quite confused about why my neutrons where reporting no energy being deposited with this method.
Thank you so much for your help and input! I’ll post some updates when I have them.
I like the find-pipe better because I can restrict the search, where grep -R will find stuff in README, .mac file, etc. Didn’t mean to offend; sorry!
If you’re using a SteppingAction or something, that’s what you have to do, since it gets called for all volumes. With an SD, the tracking essentially does that for free (the SD is only called for steps within the volume it’s attached to). G4StepStatus is an enum, so you can compare aStep->GetPreStepPoint()->GetStepStatus() with the enum value. In my own code, I typically only use it with fGeomBoundary, but there are other possibilities you can see in the .hh file.
Regarding my comment about neutrons’ energy deposit, “which is generally zero for neutrons, but might not be if you have production cuts!” In G4, the production cuts control whether secondaries get created or not. If they have energies below the production cut range, then G4 won’t actually create tracks for those secondaries, but will instead put the kinetic energy into the parent’s EnergyDeposit value.
For neutrons, this is most relevant for elastic scatters: a nucleus in a solid probably ranges out in a few tens of microns or so. Unless you set the proton (all hadrons!) production cut to zero, G4 will never create that recoiling nucleus. Instead, you’ll see an energy deposit for the neutron after the elastic scatter.
Haha no worries. It can be annoying when it digs in to Build object files, as well so find might be a better way to go.
I will check it out. Geant4/track/include/G4StepStatus.hh
This is very interesting. I hadn’t thought of a recoiling nucleus as a secondary particle. Could be quite interesting if its possible to extend G4CMP for optical phonons also using this recoil, though this is probably not necessary.
Is the scattered neutron considered a secondary?
I don’t think I will have this problem because in my sample I am only interested in how the neutrons scatter into it and their distribution within. In fact I am only interested in the depth and distribution of tritium production as a function of sample density etc. I want to record the neutron spectra and intensities at various locations where I might place the sample, in order to find the best place to do so and have some idea of the absorbed dose. Probably need to add some foils to my geometry for activation analysis, so I would have to plot the activity of the activated secondaries and take them out to a new simulation to model the gamma measurements.
Ooooh, a hidden agenda! We don’t have any plans to add optical phonons to G4CMP at this time – since downconversion scales like E^5, optical phonons decay down to a population of acoustic phonons very quickly, so we just create that population up front, by dividing a nuclear recoil Edep into Debye energy chunks. We don’t deal with defect (Frenkel) energy in G4CMP either.
Not in elastic or quasielastic scatters. For inelastic scatters (like spallation), there’s no sensible way to tag which secondary neutron “was the projectile” (darned quantum indistinguishability 
), so in G4 the neutron is killed and all of the products are secondaries.
There are some detector technologies for which G4CMP could come to be quite useful, in order to filter for signal events in the readout, apart from SuperCDMS ofcourse. Something for the future perhaps.