I’m comparing dose deposition between Geant4 and MCNP. In MCNP, it’s possible to have a volume with no material inside. For this I’m setting the density to universe_mean_density = 1.e-25 g/cm³ in Geant4. As dose deposition calculates the energy deposited per cell mass, the very low cell mass gives me weird results. I end up with the same order of magnitude of dose as in normal material cells due to the division by the tiny cell mass, despite the actual energy deposition being barely there. I have experimented a bit with other low densities, but the change in energy deposition and cell mass seem to cancel each other out.
Therefore I wanted to create a modified doseDeposit scorer which sets the dose zero if the local density is below a certain threshold. As I’m using command based scoring with a mesh, I tried to followed advice given in this thread and from a talk by M.Asai called Scoring I.
I have copied into my subroutine folder and renamed G4PSDoseDeposit, G4PSDoseDeposit3D and G4ScoreQuantityMessenger replacing “G4” with “My”. I added
G4ScoringManager *scManager = G4ScoringManager::GetScoringManager();
MyScoreQuantityMessenger *quantMess = new MyScoreQuantityMessenger(scManager);
into my main().
In MyScoreQuantityMessenger, I added a new G4UIcommand* qnoVacDoseDepCmd, copied the usual doseDeposit part in SetNewValue() and edited it in (left the cylinder mesh part as it was as I’m not going to use it). I also put the QuantityCommands() part here and again copied the usual doseDeposit stuff and replaced with my qnoVacDoseDepCmd.
When I run this, I get to the point where 1-2 threads finish (possibly without ever encountering any dose deposition since I’m trying with low history numbers) and then there’s a segmentation fault. I put some G4cout lines in and this tells me it gets to MyPSDoseDeposit(), but never to Initialize() or ProcessHits(). I haven’t even introduced the change I wanted to make to the doseDeposit scoring yet.
What am I getting wrong here? I have also tried removing all of the stuff not relating to my new scorer out of the MyScoreQuantityMessenger, but that didn’t change the results.
If you are using the command-based scorer, dose scorer will see the material of the tracking volume. This is equivalent to MCNP where you do not specify the material for the tallying volume. If you use vacuum (or ultra-thin density), there won’t be much of energy deposition.
When you use your own scorer commands, please make sure the corresponding messenger class to be instantiated for worker threads as well.
Yes, the energy deposition is low in vacuum, however, I’m comparing dose (MCNP f4 mesh tally) and the very low but non-zero density and division by cell mass causes Geant4 to have results like the upper plot (the dose deposition in the vacuum actually goes to even higher ranges which I cut off here) while MCNP produces results like the lower plot:
Therefore, I’m trying to set dose results in vacuum to zero, so that it doesn’t distract from the dose deposition in the material parts.
How do I instantiate the Messenger also for worker threads? I have tried to do it in DetectorConstruction. I need to put it into the constructor, because otherwise the new command in the input file is undefined. I have tried instantiating again in ConstructSDandField() to get it to the worker threads and then I got feedback from the MyPSDoseDeposit::ProcessHits class for the first time. Unfortunately, the segmentation fault at the end still happened. Could you explain the syntax of instantiating the Messenger for worker threads either in main() or DetectorConstruction?
ConstructSDandField() is one of the right locations to instantiate messengers for worker threads. The other location is Build() method of UserActionInitialization. Both are basically the same as long as you don’t rebuild your geometry.
Do you have the trace back to show where the code dump happened?
I tried instantiating the ScoringManager and the Messenger both in the BuildForMaster() and in the Build() part of ActionInitialization() and that seems to have done the trick. No more segmentation fault and it’s dumping results. Thank you for your help!
