| Summary: | exception with parallel world and process | ||
|---|---|---|---|
| Product: | Geant4 | Reporter: | Tom Roberts <tjrob> |
| Component: | processes/transportation | Assignee: | John Apostolakis <John.Apostolakis> |
| Status: | RESOLVED FIXED | ||
| Severity: | major | ||
| Priority: | P5 | ||
| Version: | 9.1 | ||
| Hardware: | Apple | ||
| OS: | Mac OS X | ||
| Attachments: | Patch file for G4CoupledTransportation.cc | ||
This is a more general problem than I originally thought. The exception is issued whenever G4StepLimiter limits a step right on a surface of a volume in the parallel world. The exception is actually issued at the start of the next step. While I have a workaround for the edge of the world, I have no workaround for this. And this happens much more frequently. There is a minor error in my original description: box B actually extends from z=950 to z=1050. An hour with ddd showed that the exception comes from the G4ParticleChange returned by G4CoupledTransportation::AlongStepDoit(). I also learned this additional dependency: My program has G4ElectroMagneticField::DoesFieldChangeEnergy()==true, because there can be E fields (but E is identically zero in this test). As a test I changed it to return false: the exception never happens, and all my tests succeed except those which that change invalidated. I have not gotten any test with nonzero E to work, as they all trip an unrelated bug in my code (perhaps tomorrow I'll have that fixed). Hopefully that will let you zero in on the problem more quickly, and know what your test case must include. Here's a narrative of my tracing the program: I am tracking a mu+ right down the z axis, with MaxStep=100 and with a box in the parallel world with its front face at z=100. The first step was limited by G4StepLimiter to z=100, and my UserSteppingAction printed the step. Now in G4CoupledTransportation::AlongStepGetPhysicalInteractionLength() just before the exception will be thrown, fTransportEndPosition=startPosition. Somehow startMassSafety==0.0 even though there is no boundary in the mass world anywhere nearby (there is a boundary in the parallel world right here). There is a comment: TODO: Add explicit logical status for being at a boundary, but no code. In this call which will generate the exception in the Doit(), it is right at a boundary in the parallel world. My program has DoesFieldChangeEnergy()==true, because there can be E fields. But there is none in this test. It looks to me like the underlying error is the fact that endTrackState.GetLabTimeOfFlight() returned 0, even though the particle has already moved 100 mm before this step -- hence times are negative (which caused the exception). Still in AlongStepGetPhysicalInteractionLength, endpointDistance==0.0, so it still look like this is supposed to be a 0-length step (stopped by the boundary in the parallel world, as that has not yet been crossed). safetyProposal==0 and geometryStepLength==0 also. After ASGPIL() returned its 0, the Stepping Manager was taking a 0-length step when the Doit generated the exception. After fixing the unrelated bug that prevented testing nonzero B fields (this was an error in placing rotated objects into the parallel world; downstream of a magnet I rotate coordinates and objects), I discovered that it's worse than I thought -- parallel-world-tracking does not seem to work correctly for any EM fields, electric or magnetic. When tracking through a static 0.65 Tesla B field, the tracking is just plain wrong (>100 mm difference after 2 meters, 200 MeV/c mu+, physics=QGSP). This is a comparison of my original code and my new code using the parallel world, with all objects in the parallel world artificially moved to the very end of the world so the track does not hit any of them in the region of comparison (that was to keep the steps identical for both cases, but didn't do so at all). I'm using G4EqMagElectricField for the field, and G4ClassicalRK4 for the stepper. I've been using these for several years. I still have the test version of DoesFieldChangeEnergy() returning false, and E=0 everywhere. My new code checks that there is at least one object in the parallel world, and if there are none it does not add G4ParallelWorldScoringProcess to the particles. This check happens after the parallel world is registered. When I do this, tracking is correct in the new code (i.e. the parallel world is registered, CoupledTransportation is used, but there is no G4ParallelWorldScoringProcess). Now the steps and tracking are identical (better than 0.0001 mm) for the new and the original code (still 0.65 Tesla B field). Unless I am doing something wrong here, it looks to me like the parallel world tracking is not yet ready for prime time. I don't have enough knowledge and experience of the details of Geant4 tracking to reasonably debug this further. I hope you can fix it -- tests without EM fields show this approach will speed up tracking by a factor of 10-30 (!) for simulations that use these features of my program (it currently uses an inefficient algorithm to stop at a coordinate plane, written well before parallel worlds were available). Are other people using parallel worlds with EM fields? I have been able to reproduce the problem with Geant4 9.4 patch 04, using a modified version of example/extended/field/field02, which included a Scoring Manager (triggerred with an input file) when tracking a muon.
The report is similar:
*********************************************************************************************************
* G4Track Information: Particle = e-, Track ID = 84, Parent ID = 81
*********************************************************************************************************
Step# X Y Z T KineE dEStep StepLeng TrakLeng NextVolu Process
0 -6.81 mm -38.2 cm -238 um 3.23 ns 160 keV 0 eV 0 fm 0 fm World initStep
1 -6.92 mm -39.0 cm -61.3 um 3.26 ns 80.2 MeV 1.73 keV 8.03 mm 8.03 mm World eIoni
2 -6.94 mm -40.0 cm -34.4 um 3.29 ns 179 MeV 1.62 keV 9.92 mm 1.8 cm World boxMesh_1
3 -6.94 mm -40.2 cm -31 um 3.3 ns 199 MeV 640 eV 1.93 mm 1.99 cm World eIoni
4 -6.96 mm -42.7 cm -3.12e-12 m 3.38 ns 453 MeV 3.87 keV 2.54 cm 4.53 cm World eIoni
G4ParticleChange::CheckIt : the global time goes back !!
Difference: 3.3832018586652[ns]
G4ParticleChange::CheckIt : the proper time goes back !!
Difference: 0.0038158474506449[ns]
G4ParticleChange::CheckIt
-----------------------------------------------
G4ParticleChange Information
-----------------------------------------------
# of 2ndaries : 0
-----------------------------------------------
Energy Deposit (MeV): 0
Non-ionizing Energy Deposit (MeV): 0
Track Status : Alive
True Path Length (mm) : 0
Stepping Control : 0
Mass (GeV) : 0
Charge (eplus) : 0
MagneticMoment : 0
: = 0*[e hbar]/[2 m]
Position - x (mm) : -6.96
Position - y (mm) : -427
Position - z (mm) : -3.12e-16
Time (ns) : 0
Proper Time (ns) : 0.0171
Momentum Direct - x : -0.000523
Momentum Direct - y : -1
Momentum Direct - z : 0.00129
Kinetic Energy (MeV): 453
Polarization - x : 0
Polarization - y : 0
Polarization - z : 0
Touchable (pointer) : 0x200e02a8
*** G4Exception : 200
issued by : G4ParticleChange::CheckIt
momentum, energy, and/or time was illegal
*** Event Must Be Aborted
5 -6.96 mm -42.7 cm -0.000312 fm 0 ps 453 MeV 0 eV 0 fm 4.53 cm World boxMesh_1
Created attachment 258 [details]
Patch file for G4CoupledTransportation.cc
Fixes issue with time in the test in which I reproduced the problem.
I identified a problem with the assignment of time in the case of zero steps. A patch to fix G4CoupledTransportation is attached. |
When G4StepLimiter limits a step exactly at the edge of the world volume, with a parallel world and G4ParallelWorldScoringProcess active, a fatal exception is generated. This does not happen without the parallel world process, and does not happen if CoupledTransportation limits the step at the edge of the world. So this basically only occurs when the track length from the last object to the world edge is exactly an integral multiple of the MaxStep value in the world volume, and MaxStep limits all the steps. In this example, the mass world extends to z=1350, and has Box "B" from z=1050 to z=1150; MaxStep=100 mm. The parallel world is named "ZStepWorld" and has a 1 micron thick Box with its front face at z=123.0000. The UserSteppingAction prints the name of the post-step physical volume in ZStepWorld if it is not the world PV ("ZStep PV: %s"). The G4ParallelWorldScoringProcess is named "ZStepProc". The first event has a mu+ that is 1 mrad off axis. The second event goes right down the Z axis and trips the exception. In a different run I aligned the front faces of the two physical volumes; both worlds transition into their new volumes at the same step, and inside B a 1 micron step is taken followed by a 99.9990 mm step (as expected). I have also run 100,000 events with Gaussian angles, and none of them tripped the exception; only tracks right down the Z axis do so for this example. The exception is avoided for tracks down the axis if I set MaxStep to 100.001 mm, or if I rotate B by 1 degree (CoupledTransportation then limits the step to the edge of the world, not G4StepLimiter). I do not know if the problem is in G4StepLimiter or in G4ParallelWorldScoringProcess (or possibly somewhere else, but that seems unlikely to me). Unfortunately I do not have a simple demonstration program, as I have been working within my G4beamline source (it's modularity makes that simpler than setting up a test program for things like this). ############################################################################# G4VUserPhysicsList::AddTransportation() --- G4CoupledTransportation is used ############################################################################# =================== Event 1 ================== =========== Event 1 Track 1 mu+ Parent 0 ========= Step X(mm) Y(mm) Z(mm) T(ns) KE(MeV) StepLen This Volume Process 1 0.1000 0.0000 100.0000 0.33356 120.5355 100.00000 World StepLimiter 2 0.1230 0.0000 123.0000 0.43559 120.5355 23.00006 World ZStepProc ZStep PV: Z123 3 0.1230 0.0000 123.0010 0.43559 120.5355 0.00100 World ZStepProc 4 0.2230 0.0000 223.0010 0.87917 120.5355 100.00000 World StepLimiter 5 0.3230 0.0000 323.0009 1.32274 120.5355 100.00000 World StepLimiter 6 0.4230 0.0000 423.0009 1.76632 120.5355 100.00000 World StepLimiter 7 0.5230 0.0000 523.0008 2.14357 120.5355 100.00000 World StepLimiter 8 0.6230 0.0000 623.0008 2.58714 120.5355 100.00000 World StepLimiter 9 0.7230 0.0000 723.0007 2.96439 120.5355 100.00000 World StepLimiter 10 0.8230 0.0000 823.0007 3.34165 120.5355 100.00000 World StepLimiter 11 0.9230 0.0000 923.0006 3.78522 120.5355 100.00000 World StepLimiter 12 0.9500 0.0000 950.0000 3.88708 120.5355 26.99941 World CoupledTransportation 13 1.0500 0.0000 1050.0000 4.26433 120.5355 100.00000 B StepLimiter 14 1.0500 0.0000 1050.0000 4.26433 120.5355 0.00005 B CoupledTransportation 15 1.1500 0.0000 1150.0000 4.64158 120.5355 100.00000 World StepLimiter 16 1.2500 0.0000 1249.9999 5.08515 120.5355 100.00000 World StepLimiter 17 1.3500 0.0000 1349.9999 5.52873 120.5355 100.00000 World StepLimiter 18 1.3500 0.0000 1350.0000 5.52873 120.5355 0.00015 World CoupledTransportation Event 1 Completed 1 events realTime=1 sec 1.0 ev/sec =================== Event 2 ================== =========== Event 2 Track 1 mu+ Parent 0 ========= Step X(mm) Y(mm) Z(mm) T(ns) KE(MeV) StepLen This Volume Process 1 0.0000 0.0000 100.0000 0.44358 120.5355 100.00000 World StepLimiter 2 0.0000 0.0000 123.0000 0.54560 120.5355 23.00000 World ZStepProc ZStep PV: Z123 3 0.0000 0.0000 123.0010 0.54560 120.5355 0.00100 World ZStepProc 4 0.0000 0.0000 223.0010 0.98918 120.5355 100.00000 World StepLimiter 5 0.0000 0.0000 323.0010 1.43275 120.5355 100.00000 World StepLimiter 6 0.0000 0.0000 423.0010 1.81000 120.5355 100.00000 World StepLimiter 7 0.0000 0.0000 523.0010 2.25358 120.5355 100.00000 World StepLimiter 8 0.0000 0.0000 623.0010 2.63083 120.5355 100.00000 World StepLimiter 9 0.0000 0.0000 723.0010 3.00808 120.5355 100.00000 World StepLimiter 10 0.0000 0.0000 823.0010 3.45166 120.5355 100.00000 World StepLimiter 11 0.0000 0.0000 923.0010 3.82891 120.5355 100.00000 World StepLimiter 12 0.0000 0.0000 950.0000 3.93076 120.5355 26.99900 World CoupledTransportation 13 0.0000 0.0000 1050.0000 4.30801 120.5355 100.00000 B StepLimiter 14 0.0000 0.0000 1050.0000 4.30801 120.5355 0.00000 B CoupledTransportation 15 0.0000 0.0000 1150.0000 4.68526 120.5355 100.00000 World StepLimiter 16 0.0000 0.0000 1250.0000 5.06251 120.5355 100.00000 World StepLimiter 17 0.0000 0.0000 1350.0000 5.43976 120.5355 100.00000 World StepLimiter G4ParticleChange::CheckIt : the global time goes back !! Difference: 5.439765[ns] G4ParticleChange::CheckIt : the proper time goes back !! Difference: 2.5409914[ns] G4ParticleChange::CheckIt ----------------------------------------------- G4ParticleChange Information ----------------------------------------------- # of 2ndaries : 0 ----------------------------------------------- Energy Deposit (MeV): 0 Non-ionizing Energy Deposit (MeV): 0 Track Status : Alive True Path Length (mm) : 0 Stepping Control : 0 Mass (GeV) : 2.49e-265 Charge (eplus) : 2.49e-262 MagneticMoment : 2.49e-262 : = 0*[e hbar]/[2 m] Position - x (mm) : 7.82e-316 Position - y (mm) : 0 Position - z (mm) : 1.35e+03 Time (ns) : 0 Proper Time (ns) : 0 Momentum Direct - x : -0 Momentum Direct - y : -0 Momentum Direct - z : 1 Kinetic Energy (MeV): 121 Polarization - x : 0 Polarization - y : 0 Polarization - z : 0 Touchable (pointer) : 0x9fa64c0 *** G4Exception : 200 issued by : G4ParticleChange::CheckIt momentum, energy, and/or time was illegal Severity : *** Event Must Be Aborted *** *** G4Exception: Aborting execution *** Abort trap tjrob:~/g4bl/source g4bl$