Problem 2235

Summary: G4Track returns incorrect coordinates and step lengths of chemical species
Product: Examples/Extended Reporter: Ivan <beo0>
Component: medical/dnaAssignee: hoangtran <ngoc.hoang.tran>
Status: RESOLVED FIXED    
Severity: normal CC: ngoc.hoang.tran
Priority: P4    
Version: 10.6   
Hardware: PC   
OS: Windows   
Attachments: TimeStepAction.cc and TimeStepAction.hh files
Macro file to use with "chem6" example of Geant4 version 10.7-beta

Description Ivan 2020-04-07 19:36:53 CEST
Created attachment 613 [details]
TimeStepAction.cc and TimeStepAction.hh files

Hello

I would like to report a problem with DNA chemistry module. I was trying to model paths of all chemical secies in water and found two problems that could be bugs. First one is minor: global times related to chemical species at the same time point can slightly differ. This difference is very small and looks like it's just an error related to a limited accuracy of storing numbers as "double" type in C++. However, times are formaly not equal. I don't know if this is important, so just saying... Second problem is more serious. The problem is that some chemical species are stuck between "pre" time step and "post" time step. So they do not move, thier global times does not change as well as thier coordinates. But at the same time the function track->GetStepLength() returns non-zero step telling user that specie actually moves between pre time step and post time step. It also happens that specie is not frozen, but GetStepLength() still gives value that is sufficiently different from the value calculated by coordinates of specie at pre step and post step. Such mismatches occur pretty rare, probably in 1% of time steps or even more rare. Thus the effect from those probably cannot be seen in general results produced by DNA chemistry, e.g. when user calculates total number of species. However, the impact of this on results which take into account details of species tracks is unclear.

How to reproduce this bugs:
Take "chem2" example and replace TimeStepAction.cc and TimeStepAction.hh files with ones I attached to this report. To obtain more readable output it might be a good idea to model other primary particle: replace 50 keV proton with 1 keV electron in beam.in file. Run resulting example and make it output everything into text file instead of console. For example, on Windows PowerShell one can write command to execute example (from build directory):

Debug\./chem2 > tmp.txt

Then open tmp.txt file and search for "WARNING" or "ERROR" messages (ctrl+f).

The program stores track IDs, coordinates, kinetic energies, molecule IDs and global times of all chemical species at pre time step and post time step and checks consistency of stored data. It checks:
1) if properties of tracks with the same ID were changed between post step of previous step and pre step of current step
2) if tracks were added or removed between post step of previous step and pre step of current step
3) if molecule IDs of tracks with the same IDs match at pre step and post step
4) if global times of species are the same both for pre step and post step
5) if step length returned by GetStepLength() matches with step length calculated from coordinates of chemical specie
6) if step length equals zero

For now I performed only few tests and I didn't observe any problems in 1,2 and 3, but 4,5 and 6 had problems.

I would appreciate any info of how to fix the problem or if it can be neglected.

Best Regards,
Ivan
Comment 1 hoangtran 2020-04-20 08:18:51 CEST
the bug is being investigated by Mathieu and Hoang. I guest the problem may come from the step to geometry boundary. However, in chemistry, we care more timestep than StepLength (because it's brownian motion) so if you want to know about position please use GetPosition(). 

Best regards,
Hoang
Comment 2 Sebastien Incerti 2020-07-16 23:31:29 CEST
Reassigned
Comment 3 Ivan 2020-07-23 13:31:25 CEST
Created attachment 635 [details]
Macro file to use with "chem6" example of Geant4 version 10.7-beta
Comment 4 Ivan 2020-07-23 13:33:06 CEST
Comment on attachment 635 [details]
Macro file to use with "chem6" example of Geant4 version 10.7-beta

Hello

I would like to report another problem that occurs in "chem" examples. It might be related to the same bug (see the first post), so I've decided to not create new topic. Unlike previous bug this one completely eliminates the possibility of getting results of modeling, because program terminates with an error(G4Exception). Exception occurs in both 10.6 and 10.7-beta versions of Geant4.

How to reproduce: take "chem5" or "chem6" example of 10.6 or 10.7-beta version of Geant4 and change size of "world volume" from 1 km to 5 um. Then open "beam.in" file and change source to alpha particles with energy of 40.49 MeV and remove the possibility of track being killed by "primaryKiller". For "chem6" example you can use "beam.in" file attached to this message. Note that other types of particles also produce this exception, but less frequently.

Unfortunately, it might take several days of calculation before the exception occurs, but for alpha particles it is sure to occur.

The exception occurs for the following constructors: G4EmDNAChemistry, G4EmDNAChemistry_option1 and G4EmDNAChemistry_option2. It does not occur for G4EmDNAChemistry_option3, which is very fast, but doesn't allow user to access information of tracks of chemical species on each time step (using G4Track class at pre step and post step in TimeStepAction.cc file, see the first message). It also gives G-values, that might be significally different from ones when using other chemistry constructors.

This is how exception looks like:


-------- EEEE ------- G4Exception-START -------- EEEE -------
*** G4Exception : ITManager002
      issued by : G4ITManager::FindNearest
Bad request : no node found in the IT you are searching closest neighbourg for
*** Fatal Error In Argument *** core dump ***
G4Track (00000258C2DA7090) - track ID = -10882, parent ID = -10530
 Particle type : H_2 - creator process : H2O_DNAMolecularDecay, creator model : Undefined
 Kinetic energy : 0.0385195 eV  - Momentum direction : (-0.749714,-0.340201,-0.56762)
 Step length : 1.52532e-297 fm  - total energy deposit : 7.35921 eV
 Pre-step point : (-5.68754e-08,-1.37734e-06,-0.00249829) - Physical volume : World (G4_WATER)
 - defined by : e-_G4DNAVibExcitation - step status : 4
 Post-step point : (-5.68754e-08,-1.37734e-06,-0.00249829) - Physical volume : World (G4_WATER)
 - defined by : e-_G4DNAElectronSolvation - step status : 4
 *** Note: Step information might not be properly updated.

-------- EEEE -------- G4Exception-END --------- EEEE -------


*** G4Exception: Aborting execution ***


Coordinates of chemical specie in the exception are always close to boundary of the world volume. I tried to use track->SetTrackStatus(fStopAndKill) to kill all species, that are close to boundary, but it resulted in another exception, so this doesn't solve the problem.
Comment 5 hoangtran 2021-08-25 17:18:55 CEST
Dear Ivan, 

Thank you very much for your report! 

Sorry for relying late!

Please don't change the size of "world volume" to 5 um. The Brownian motion close to boundary may have problem. We are working on this issue. 

For G4EmDNAChemistry_option3, it used another approach called IRT which doesn't track the molecules so we don't have information of tracks to access.


Best regards,
Hoang
Comment 6 hoangtran 2022-10-07 10:44:46 CEST
This issue should be fixed on G4.11.1.