Geant4 Version: 11.2.1
Operating System: Red Hat Enterprise Linux 9.4
Compiler/Version: GNU 12.2
In my simulation of a capture beamline after a thick converter, I regularly observe code crashes when the particles propagate through a magnetic field. I can reproduce these errors with the following in my void DetectorConstruction::SetupGeometry() :
G4double converterlength = 5*mm;
auto world_box = new G4Box("world_box", 2.5*m, 2.5*m, 5.*m);
auto world_log = new G4LogicalVolume(world_box, exphall_material, "world_log", 0, 0, 0, false);
new G4PVPlacement(0, G4ThreeVector(), world_log, "World",0, false, 0);
auto converter_box = new G4Box("converter_box",.04*m,.04*m, converterlength);
auto converter_log = new G4LogicalVolume(converter_box, converter_material, "converter_log",0,0,0,false);
new G4PVPlacement(0, G4ThreeVector(0.*m, 0.*m, -converterlength),
converter_log, "converter", world_log, false, 0);
// Make the Bfield volume and field
auto B_box = new G4Box("Bfield", 5*cm, 5*cm, 10*cm);
auto B_log = new G4LogicalVolume(B_box, exphall_material, "Bfield_log", 0, 0, 0);
new G4PVPlacement(0, G4ThreeVector(0., 0., 15*cm), B_log, "Bfield_phys", world_log, 0, 0);
auto BField = new G4QuadrupoleMagField(100*tesla/m);
auto localFieldManager = new G4FieldManager();
localFieldManager->SetDetectorField(BField);
auto pEquation = new G4Mag_UsualEqRhs(BField);
auto pStepper = new G4DormandPrince745(pEquation);
G4double stepMinimum = .01 * millimeter;
auto pChordFinder = new G4ChordFinder(BField, stepMinimum, pStepper);
localFieldManager->SetChordFinder( pChordFinder );
B_log->SetFieldManager(localFieldManager, true);
where exphall_material is G4_AIR and the converter is G4_W. The physics list is made with auto physicsList = physListFactory->GetReferencePhysList("FTFP_BERT");.
When running the code, I see two types of issues:
- An exception/warning
-------- WWWW ------- G4Exception-START -------- WWWW -------
*** G4Exception : GeomField1001
issued by : G4IntegrationDriver::AccurateAdvance()
Proposed step is zero; hstep = 0 !
*** This is just a warning message. ***
-------- WWWW -------- G4Exception-END --------- WWWW -------
which seems to be similar to what was reported for Belle simulations by Dennis Wright a few years ago. These are annoying, but seemingly harmless (?).
- An exception, leading to a core dump:
G4WT93 >
-------- EEEE ------- G4Exception-START -------- EEEE -------
*** G4Exception : GeomNav0003
issued by : G4MultiLevelLocator::EstimateIntersectionPoint()
Current Position and Direction
Step# s X(mm) Y(mm) Z(mm) N_x N_y N_z Delta|N| StepLen StartSafety PhsStep
Step taken was -2.26763e-10 out of PhysicalStep= -1
Final safety is: 0.0065435
Chord length = 0.00836269
Error in advancing propagation.
The final curve point is NOT further along than the original!
Going *backwards* from len(A) = 65.3974 to len(B) = 65.3974
Curve distance is -2.26763e-10 mm
Point A' (start) is ( X= 13.1191581 -13.9091056 50.0065435 P= 10.1289694 -2.11995412 13.4693282 Pmag= 16.9856713 Ekin= 16.4823571 l= 65.3974084809 m0= 0.510999 (Pdir-1)= 0 t_lab= 0 t_proper= 0 PolV= (0,0,0) )
Point B' (end) is ( X= 13.114488 -13.9076678 49.9997569 P= 10.1319531 -2.11662265 13.4671655 Pmag= 16.9853205 Ekin= 16.4820065 l= 65.3974084807 m0= 0.510999 (Pdir-1)= -1.11022e-16 t_lab= 0 t_proper= 0 PolV= (0,0,0) )
fEpsStep= 6.6863e-05
In full precision, the position, momentum, E_kin, length, rest mass ... are:
Point A[0] (Curve start) is ( X= 0.858643058 16.8575348 53.9024085 P= 0.796854113 8.55022088 35.7380418 Pmag= 36.7552566 Ekin= 36.2478096 l= 64.5730527357 m0= 0.510999 (Pdir-1)= 0 t_lab= 0 t_proper= 0 PolV= (0,0,0) )
Point S (Sub start) is ( X= 18.0729387 45.0931841 249.98751 P= 46.2533607 -6.60168169 108.406237 Pmag= 118.046041 Ekin= 117.536148 l= 65.3953799373 m0= 0.510999 (Pdir-1)= 0 t_lab= 0 t_proper= 0 PolV= (0,0,0) ) Point A' (Current start) is ( X= 13.1191581 -13.9091056 50.0065435 P= 10.1289694 -2.11995412 13.4693282 Pmag= 16.9856713 Ekin= 16.4823571 l= 65.3974084809 m0= 0.510999 (Pdir-1)= 0 t_lab= 0 t_proper= 0 PolV= (0,0,0) )
Point E (Trial Point) is (13.114655264755576525,-13.907719265967845246,49.999999999999992895)
Point F (Intersection) is ( X= 13.1191581 -13.9091056 50.0065435 P= 10.1289694 -2.11995412 13.4693282 Pmag= 16.9856713 Ekin= 16.4823571 l= 65.3974084809 m0= 0.510999 (Pdir-1)= 0 t_lab= 0 t_proper= 0 PolV= (0,0,0) )
Point B' (Current end) is ( X= 13.114488 -13.9076678 49.9997569 P= 10.1319531 -2.11662265 13.4671655 Pmag= 16.9853205 Ekin= 16.4820065 l= 65.3974084807 m0= 0.510999 (Pdir-1)= -1.11022e-16 t_lab= 0 t_proper= 0 PolV= (0,0,0) )
Point B[0] (Curve end) is ( X= 18.1148767 45.0819149 250.045438 P= 46.3883077 -6.94116896 108.328935 Pmag= 118.047503 Ekin= 117.53761 l= 66.7277570832 m0= 0.510999 (Pdir-1)= 0 t_lab= 0 t_proper= 0 PolV= (0,0,0) )
LocateIntersection parameters are :
Substep no (total) = 3
Substep no = 1 at depth= 1
* Location: G4MultiLevelLocator::EstimateIntersectionPoint()- After EndIf(Intersects_AF)
* Bool flags: Recalculated = 0 Intersects_AF = 0 Intersects_FB = 1
* Number of calls to MLL:EIP= 1893
*** Fatal Exception *** core dump ***
G4WT93 > G4Track (0x149af0301e80) - track ID = 106, parent ID = 26
G4WT93 > Particle type : e- - creator process : conv, creator model : model_EplusEminisPair
G4WT93 > Kinetic energy : 61.2916 MeV - Momentum direction : (0.398306,-0.0596073,0.915314)
G4WT93 > Step length : 5.56865 cm - total energy deposit : 0 eV
G4WT93 > Pre-step point : (17.3462,44.4346,248.027) - Physical volume : Bfield_phys (G4_AIR)
G4WT93 > - defined by : eIoni - step status : 4
G4WT93 > Post-step point : (17.3462,44.4346,248.027) - Physical volume : Bfield_phys (G4_AIR)
G4WT93 > - defined by : eIoni - step status : 7
G4WT93 > *** Note: Step information might not be properly updated.
G4WT93 >
-------- EEEE -------- G4Exception-END --------- EEEE -------
This happens on different threads at different times, but never allows simulating the setup properly through. It most often happens at the boundary of the field, but sometimes also inside the field. I’ve tried changing the DeltaChord, StepMinimum, apply G4UserLimits, try different steppers, change eps etc, but this core dump always ends up happening. Sometimes it’s relativistic electrons/positrons (10-100MeV), sometimes low energy electrons/positrons, but it always results in a crash.
Interestingly, same code, same G4 version, but running on MacOS with M1 and CLang 15.0.0 never ends up crashing: the same exception is shown, but the code carries on.
Does anyone have any idea what this could be caused by or where I should look to try to make it better?
Best,
Kris