When using the G4PhantomParameterisation, a bad_alloc occurs when a particle moves very close and parallel to a voxel edge. The bad_alloc occures in the for( ; ; ) loop in G4RegularNavigation::ComputeStepSkippingEqualMaterials, because DistanceToOut( localPoint, localDirection ) of G4Box returns 0 and therefore the length is only incremented by kCarTolerance resulting in many loop iterations. In each iteration, G4RegularNavigationHelper::Instance()->AddStepLength( copyNo, newStep ) is called, adding an element to a std::vector until the bad_alloc occurs. This can be reproduced usig e.g. localPosition = (-0.5,-0.154576,0.5), localDirection = (1,0,6.12323e-17) and a voxel halfsize of (0.5,0.5,0.5) I see, that this problem occurs due to floating point precision issues and the root cause (G4Box returing 0) might be hard to fix without maybe breaking assumptions elsewhere. Therefore I propose to change the function G4RegularNavigation::ComputeStepSkippingEqualMaterials such, that an additional check is added, to see if the copyNo has changed (which is assumed implicitly at the moment) and sum up all lengths in the same copyNo before adding it to the G4RegularNavigationHelper::theStepLengths vector. This would result in a single slow event but prevent the memory overflow
I was not able to reproduce your problem with the DICOM examples and using a simple test.g4dcm (see below). As you say the precision is a delicate point in G4RegularNavigation.cc And I do not understand how your solution solves the problem: it will not create more StepLength's, but it will be kept looping forever. A safer solution would be to push the step a minimal amount when it is 0, but I would like to be able to reproduce your problem to test it. test.g4dcm =========== 1 0 Water 10 10 10 -5 5 -5 5 -5 5 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.
For what its worth, it already pushes the step a minimal amount when voxelBox->DistanceToOut returns zero: line 209 of G4RegularNavigation pushes it an additional kCarTolerance. But if the track is nearly parallel to the voxel boundary, it could take a very large number of such steps before it either is far enough from the boundary that DistanceToOut doesn't return zero, or the for loop breaks because it enters a new material or it leaves the parameterized volume entirely. I think one potential workaround until this bug is fixed is to call parameterization->SetSkipEqualMaterials(false); this will cause the navigator to ultimately use G4NormalNavigation::ComputeStep to calculate the step. This way, when a track is close to & parallel to a voxel boundary, G4Navigator gets a chance to check for 0 length steps and can start moving it 100*kCarTolerance if it appears stuck, eventually killing the track if that doesn't work. In the case of the bug, G4Navigator doesn't get the chance to move or kill such tracks since ComputeStepSkippingEqualMaterials doesn't return a step size to G4Navigator until the for (;;) loop breaks (and since the step that it returns will be at minimum kCarTolerance, it wouldn't trigger the zero length step code in G4Navigator anyways, as it looks for steps of size < 0.05*kCarTolerance). Thus, the solution would need to be logic inside that for loop. I believe koernigchris's solution should work: tracks that are parallel to and very close to a voxel boundary would run slowly (potentially moving through the entire parameterised volume one kCarTolerance at at time if the proposed step length is large enough and the track leads through all the same material type), but it would only generate a new vector entry each time they enter a new copyNo so it wouldn't run out of memory in the process of all these iterations. Another option would be to implement logic similar to G4Navigator re: multiple zero length steps, first trying to move more that kCarTolerance at a time, then simply killing the track if it still remains stuck.
(In reply to Pedro.Arce from comment #1) > I was not able to reproduce your problem with the DICOM examples and using a > simple test.g4dcm (see below). As you say the precision is a delicate point > in G4RegularNavigation.cc > And I do not understand how your solution solves the problem: it will not > create more StepLength's, but it will be kept looping forever. > A safer solution would be to push the step a minimal amount when it is 0, > but I would like to be able to reproduce your problem to test it. > > > test.g4dcm > =========== > 1 > 0 Water > 10 10 10 > -5 5 > -5 5 > -5 5 > 0 0 0 0 0 0 0 0 0 0 > 0 0 0 0 0 0 0 0 0 0 > 0 0 0 0 0 0 0 0 0 0 > 0 0 0 0 0 0 0 0 0 0 > 0 0 0 0 0 0 0 0 0 0 > 0 0 0 0 0 0 0 0 0 0 > 0 0 0 0 0 0 0 0 0 0 > 0 0 0 0 0 0 0 0 0 0 > 0 0 0 0 0 0 0 0 0 0 > 0 0 0 0 0 0 0 0 0 0 > 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. > 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. > 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. > 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. > 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. > 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. > 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. > 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. > 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. > 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. Alright, I think I have all the steps you'll need to reproduce this bug using the DICOM example and your test.g4dcm: * First line of your Data.dat should be 1 (this way, with your g4dcm file as it is, you will have 1mmx1mmx1mm voxels) * In DicomPhantomParameterizationColour.cc line 53: SetSkipEqualMaterials(true); //(this bug occurs in G4RegularNavigation::ComputeStepSkippingEqualMaterials, but the DICOM example normally sets this parameter to false, so it won't actually enter the bugged code) * In DicomPrimaryGeneratorAction.cc line 81: G4ThreeVector dir(1e-8,0,1); //Very close to parallel, with a small positive x component so that p.x()*v.x() > 0 * In the same file line 87: fParticleGun->SetParticlePosition(G4ThreeVector(0.99999999999*mm,0.,0)); //Very close to the wall of the voxel, positioned so that with the small positive x component p.x()*v.x() > 0 Run it, and when the default vis.mac file runs the beamOn, memory usage will ramp up very quickly. On my machine with the example otherwise left to the default, it'll eventually throw a bad_alloc once the process memory approaches around 99GB, this probably depends on the number of threads and the amount of memory available on the particular machine used.
Created attachment 596 [details] G4RegularNavigation.cc
Created attachment 597 [details] G4RegularNavigation.hh
I worte a solution as you suggested, and in the style of G4Navigator. Please test if it works. And thanks
Created attachment 637 [details] G4RegularNavigation.cc
Sorry for my very late reply and the double post - I am not familiar with bugzilla Unfortunatly, this solution does not work for me. It still get momory overflow and it is very verbose. I think this post on the forum is also related to this issue: https://geant4-forum.web.cern.ch/t/g4phantomparameterisation-and-stuck-particles/2275 I have an alternative solution, which is based on the previous version of the code and works well for me. I also found, that the DICOM example can be used for testing, when reducing the voxel size to around 0.1mm and adding some spread in the x-y position of the generated particle. I also changed the particle to gamma at 100keV to match my use case but i guess it should also be reproducable with electrons.
