How would you simulate an irregularly shaped voxel grid, like the red volume in the figure below:
My current implementation consists of a G4Box as the “root” volume, using a G4VNestedParameterisation to assign different materials to the voxels. Empty voxels are assigned the background material (ie air). For scoring, I use a custom class derived from G4VSensitiveDetector.
Now the issue is that I would like that the space of the empty voxels could be occupied by some other volume. I can think of the following options:
Subtracting a chunk from the G4Box root solid.
The documentation says that replicas must completely fill their containing mother volume but it does not say that some of the replicas can be outside. I did a quick check, replacing the root G4Box with a G4Tubs, and it seems to work and pass the geometry checks.
Option 1: subtracting the blue solid. The results which were scored in voxels that were chopped would need weighing to make sense, but luckily I don’t need them.
Is there any simple way to check that only empty voxels are affected?
Option 2: subtracting the empty voxels. This seems like a bad idea, unless something inordinately smart is done: the number of empty voxels can be large, and there are lot of matching surfaces. Furthermore, this may end up creating disjoined solids, which is not a valid Boolean composition (side question: what does this mean? Crash, meaningless results?).
Give up on the current approach, and create one physical volume and one detector for each non-empty voxel.
The amount of detectors would be huge, and the per-thread memory footprint would be prohibitive.
Is there any other possibility that I am missing? Any other piece of wisdom you would want to share?
Ah, of course, that would be a possibility for the detector!
There’s an issue I observed with parallel worlds which I posted a while back: in summary, my simulations misbehaved if the parallel world contained a volume right on top of an identical volume in the physical world. I did not investigate the issue further (I was lazy and added a small displacement…), but it looked as if the parallel navigator was getting confused when two surfaces were crossed at the same instant. If this is the case, that would limit this approach.
Give a try deleting the SkipEqualMaterial option of parallel navigation.
Alternatively you may use partialDICOM (see DICOM example), although in my experience this is very slow
I am not using G4PhantomParameterisation, but my own parameterisation derived from G4VNestedParameterisation.
I am trying to understand the parameterisations used in the DICOM example. Am I understanding correctly that their aim is efficient navigation? Would they address the overlapping problem I discussed? As far as I can tell, the container volume is a G4Tubs constructed in DicomPartialDetectorConstruction::ConstructPhantomContainer(). If this volume overlaps with some other volume, this would cause undefined behavior, wouldn’t it?
The G4Tubs is just an example. You can use the volume you want to contain the partial DICOM. Just be sure the .g4pdcm file voxels do not protrude from it. And, as you say the container volume does not overlap with any other
Hm, but that is my problem: in my geometric setup, I need to allow for some volume overlapping with the empty voxels. Does the partial DICOM approach support this?
With partial DICOM the empty volume do not exist. if the track exits the non-empty volumes it exits the DICOM. The only way to allow for volume overlapping in Geant4 is with the parallel worlds.
I am trying to understand how the partial DICOM achieves it. Is it because G4PartialPhantomParameterisation::ComputeVoxelIndices() will only map copy numbers to filled voxel coordinates?
Does this mean that with G4PVParameterised daughter volumes are not required to fill the mother volume exactly, and they must not protrude outside of it?
@arce I tested a voxel implementation using G4PVParameterised (similar to G4PartialPhantomParameterisation), and you are right: it is woefully slow, and not really a viable option.
With G4PVParameterised you have to be careful with the axis, kXaxis instead of kUndefined. With the bad option it is very slow because the Geant4 navigation SmartVoxel optimization is needed to make it fast (Geant4 internally divides the world in voxels so that when a track exits a volume it only loops for parents amont the volumes in the next voxel)
Thank you so much for your patience… I am finally starting to understand how it works. So, if I want to use the efficient navigation of G4RegularNavigator I should use a parameterisation derived from G4PhantomParameterisation, and call SetRegularStructureId(i) with i>0 for the parameterised voxel physical volume.
By the way, does the value passed to SetRegularStructureId matter, other than being larger than 0 or not? Is it left like that for future changes, or how is it used in user applications?
I recognise I do not recall why it is like that, but it must be so, because in
G4GeometryManager.cc
&& (volume->GetDaughter(0)->GetRegularStructureId()!=1) ) )
It is a code I wrote more than 10 years ago, I would some time to remember this.
In case someone ends up landing here: the way to implement these overlapping volumes is using layered mass geometry. Probably @allison was suggesting this, but I blindly ignored that section in the documentation, since I was wrongly convinced that parallel world had to be massless.
Besides the documentation, there is this presentation by @makotoasai, where radiotherapy applications are explicitly mentioned.
There is also an article which might be of interest.
