Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/electromagnetic/dna/models/src/G4DNAMoleculeEncounterStepper.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

Diff markup

Differences between /processes/electromagnetic/dna/models/src/G4DNAMoleculeEncounterStepper.cc (Version 11.3.0) and /processes/electromagnetic/dna/models/src/G4DNAMoleculeEncounterStepper.cc (Version 11.1)


  1 //                                                  1 //
  2 // *******************************************      2 // ********************************************************************
  3 // * License and Disclaimer                         3 // * License and Disclaimer                                           *
  4 // *                                                4 // *                                                                  *
  5 // * The  Geant4 software  is  copyright of th      5 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
  6 // * the Geant4 Collaboration.  It is provided      6 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
  7 // * conditions of the Geant4 Software License      7 // * conditions of the Geant4 Software License,  included in the file *
  8 // * LICENSE and available at  http://cern.ch/      8 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
  9 // * include a list of copyright holders.           9 // * include a list of copyright holders.                             *
 10 // *                                               10 // *                                                                  *
 11 // * Neither the authors of this software syst     11 // * Neither the authors of this software system, nor their employing *
 12 // * institutes,nor the agencies providing fin     12 // * institutes,nor the agencies providing financial support for this *
 13 // * work  make  any representation or  warran     13 // * work  make  any representation or  warranty, express or implied, *
 14 // * regarding  this  software system or assum     14 // * regarding  this  software system or assume any liability for its *
 15 // * use.  Please see the license in the file      15 // * use.  Please see the license in the file  LICENSE  and URL above *
 16 // * for the full disclaimer and the limitatio     16 // * for the full disclaimer and the limitation of liability.         *
 17 // *                                               17 // *                                                                  *
 18 // * This  code  implementation is the result      18 // * This  code  implementation is the result of  the  scientific and *
 19 // * technical work of the GEANT4 collaboratio     19 // * technical work of the GEANT4 collaboration.                      *
 20 // * By using,  copying,  modifying or  distri     20 // * By using,  copying,  modifying or  distributing the software (or *
 21 // * any work based  on the software)  you  ag     21 // * any work based  on the software)  you  agree  to acknowledge its *
 22 // * use  in  resulting  scientific  publicati     22 // * use  in  resulting  scientific  publications,  and indicate your *
 23 // * acceptance of all terms of the Geant4 Sof     23 // * acceptance of all terms of the Geant4 Software license.          *
 24 // *******************************************     24 // ********************************************************************
 25 //                                                 25 //
 26 //                                                 26 //
 27 // Author: Mathieu Karamitros (kara@cenbg.in2p     27 // Author: Mathieu Karamitros (kara@cenbg.in2p3.fr)
 28 //                                                 28 //
 29 // WARNING : This class is released as a proto     29 // WARNING : This class is released as a prototype.
 30 // It might strongly evolve or even disapear i     30 // It might strongly evolve or even disapear in the next releases.
 31 //                                                 31 //
 32 // History:                                        32 // History:
 33 // -----------                                     33 // -----------
 34 // 10 Oct 2011 M.Karamitros created                34 // 10 Oct 2011 M.Karamitros created
 35 //                                                 35 //
 36 // -------------------------------------------     36 // -------------------------------------------------------------------
 37                                                    37 
 38 #include "G4DNAMoleculeEncounterStepper.hh"        38 #include "G4DNAMoleculeEncounterStepper.hh"
 39                                                <<  39 #include "G4VDNAReactionModel.hh"
 40 #include "G4DNAMolecularReactionTable.hh"          40 #include "G4DNAMolecularReactionTable.hh"
 41 #include "G4H2O.hh"                                41 #include "G4H2O.hh"
 42 #include "G4MolecularConfiguration.hh"         << 
 43 #include "G4MoleculeFinder.hh"                 << 
 44 #include "G4UnitsTable.hh"                     << 
 45 #include "G4VDNAReactionModel.hh"              << 
 46 #include "G4memory.hh"                             42 #include "G4memory.hh"
 47                                                <<  43 #include "G4UnitsTable.hh"
 48 #include <memory>                              <<  44 #include "G4MoleculeFinder.hh"
                                                   >>  45 #include "G4MolecularConfiguration.hh"
 49                                                    46 
 50 using namespace std;                               47 using namespace std;
 51 using namespace CLHEP;                             48 using namespace CLHEP;
 52                                                    49 
 53 //#define DEBUG_MEM                                50 //#define DEBUG_MEM
 54                                                    51 
 55 #ifdef DEBUG_MEM                                   52 #ifdef DEBUG_MEM
 56 #include "G4MemStat.hh"                            53 #include "G4MemStat.hh"
 57 using namespace G4MemStat;                         54 using namespace G4MemStat;
 58 #endif                                             55 #endif
 59                                                    56 
 60 G4DNAMoleculeEncounterStepper::Utils::Utils(co     57 G4DNAMoleculeEncounterStepper::Utils::Utils(const G4Track& tA,
 61                                             co     58                                             const G4MolecularConfiguration* pMoleculeB)
 62     : fpTrackA(tA)                                 59     : fpTrackA(tA)
 63     , fpMoleculeB(pMoleculeB)                      60     , fpMoleculeB(pMoleculeB)
 64 {                                                  61 {
 65     fpMoleculeA = GetMolecule(tA);                 62     fpMoleculeA = GetMolecule(tA);
 66     fDA = fpMoleculeA->GetDiffusionCoefficient     63     fDA = fpMoleculeA->GetDiffusionCoefficient();
 67     fDB = fpMoleculeB->GetDiffusionCoefficient     64     fDB = fpMoleculeB->GetDiffusionCoefficient();
 68     fConstant = 8 * (fDA + fDB + 2 * sqrt(fDA      65     fConstant = 8 * (fDA + fDB + 2 * sqrt(fDA * fDB));
 69 }                                                  66 }
 70                                                    67 
 71 G4DNAMoleculeEncounterStepper::G4DNAMoleculeEn     68 G4DNAMoleculeEncounterStepper::G4DNAMoleculeEncounterStepper()
 72     :                                          <<  69     : G4VITTimeStepComputer()
 73      fMolecularReactionTable(reference_cast<co <<  70     , fHasAlreadyReachedNullTime(false)
                                                   >>  71     , fMolecularReactionTable(reference_cast<const G4DNAMolecularReactionTable*>(fpReactionTable))
                                                   >>  72     , fReactionModel(nullptr)
                                                   >>  73     , fVerbose(0)
 74 {                                                  74 {
 75     fpTrackContainer = G4ITTrackHolder::Instan     75     fpTrackContainer = G4ITTrackHolder::Instance();
 76     fReactionSet = G4ITReactionSet::Instance()     76     fReactionSet = G4ITReactionSet::Instance();
 77 }                                                  77 }
 78                                                    78 
 79 G4DNAMoleculeEncounterStepper::~G4DNAMoleculeE     79 G4DNAMoleculeEncounterStepper::~G4DNAMoleculeEncounterStepper() = default;
 80                                                    80 
 81 void G4DNAMoleculeEncounterStepper::Prepare()      81 void G4DNAMoleculeEncounterStepper::Prepare()
 82 {                                                  82 {
 83     G4VITTimeStepComputer::Prepare();              83     G4VITTimeStepComputer::Prepare();
 84                                                    84 
 85 #if defined (DEBUG_MEM)                            85 #if defined (DEBUG_MEM)
 86     MemStat mem_first, mem_second, mem_diff;       86     MemStat mem_first, mem_second, mem_diff;
 87 #endif                                             87 #endif
 88                                                    88 
 89 #if defined (DEBUG_MEM)                            89 #if defined (DEBUG_MEM)
 90     mem_first = MemoryUsage();                     90     mem_first = MemoryUsage();
 91 #endif                                             91 #endif
 92     G4MoleculeFinder::Instance()->UpdatePositi     92     G4MoleculeFinder::Instance()->UpdatePositionMap();
 93                                                    93 
 94 #if defined (DEBUG_MEM)                            94 #if defined (DEBUG_MEM)
 95     mem_second = MemoryUsage();                    95     mem_second = MemoryUsage();
 96     mem_diff = mem_second - mem_first;             96     mem_diff = mem_second - mem_first;
 97     G4cout << "\t || MEM || G4DNAMoleculeEncou     97     G4cout << "\t || MEM || G4DNAMoleculeEncounterStepper::Prepare || "
 98         "After computing G4ITManager<G4Molecul     98         "After computing G4ITManager<G4Molecule>::Instance()->"
 99         "UpdatePositionMap, diff is : " << mem     99         "UpdatePositionMap, diff is : " << mem_diff << G4endl;
100 #endif                                            100 #endif
101 }                                                 101 }
102                                                   102 
103 void G4DNAMoleculeEncounterStepper::Initialize    103 void G4DNAMoleculeEncounterStepper::InitializeForNewTrack()
104 {                                                 104 {
105     if (fReactants)                               105     if (fReactants)
106     {                                             106     {
107         fReactants.reset();                       107         fReactants.reset();
108     }                                             108     }
109     fSampledMinTimeStep = DBL_MAX;                109     fSampledMinTimeStep = DBL_MAX;
110     fHasAlreadyReachedNullTime = false;           110     fHasAlreadyReachedNullTime = false;
111 }                                                 111 }
112                                                   112 
113 template<typename T>                              113 template<typename T>
114 inline bool IsInf(T value)                        114 inline bool IsInf(T value)
115 {                                                 115 {
116     return std::numeric_limits<T>::has_infinit    116     return std::numeric_limits<T>::has_infinity
117         && value == std::numeric_limits<T>::in    117         && value == std::numeric_limits<T>::infinity();
118 }                                                 118 }
119                                                   119 
120 G4double                                          120 G4double
121 G4DNAMoleculeEncounterStepper::CalculateStep(c    121 G4DNAMoleculeEncounterStepper::CalculateStep(const G4Track& trackA,
122                                              c    122                                              const G4double& userMinTimeStep)
123 {                                                 123 {
124     auto pMoleculeA = GetMolecule(trackA);        124     auto pMoleculeA = GetMolecule(trackA);
125     InitializeForNewTrack();                      125     InitializeForNewTrack();
126     fUserMinTimeStep = userMinTimeStep;           126     fUserMinTimeStep = userMinTimeStep;
127                                                   127 
128 #ifdef G4VERBOSE                                  128 #ifdef G4VERBOSE
129     if (fVerbose != 0)                         << 129     if (fVerbose)
130     {                                             130     {
131         G4cout                                    131         G4cout
132             << "______________________________    132             << "_______________________________________________________________________"
133             << G4endl;                            133             << G4endl;
134         G4cout << "G4DNAMoleculeEncounterStepp    134         G4cout << "G4DNAMoleculeEncounterStepper::CalculateStep" << G4endl;
135         G4cout << "Check done for molecule : "    135         G4cout << "Check done for molecule : " << pMoleculeA->GetName()
136             << " (" << trackA.GetTrackID() <<     136             << " (" << trackA.GetTrackID() << ") "
137             << G4endl;                            137             << G4endl;
138     }                                             138     }
139 #endif                                            139 #endif
140                                                   140 
141     //________________________________________    141     //__________________________________________________________________
142     // Retrieve general informations for makin    142     // Retrieve general informations for making reactions
143     auto pMolConfA = pMoleculeA->GetMolecularC    143     auto pMolConfA = pMoleculeA->GetMolecularConfiguration();
144                                                   144 
145     const auto pReactantList = fMolecularReact    145     const auto pReactantList = fMolecularReactionTable->CanReactWith(pMolConfA);
146                                                   146 
147     if (pReactantList == nullptr)              << 147     if (!pReactantList)
148     {                                             148     {
149 #ifdef G4VERBOSE                                  149 #ifdef G4VERBOSE
150         //    DEBUG                               150         //    DEBUG
151         if (fVerbose > 1)                         151         if (fVerbose > 1)
152         {                                         152         {
153             G4cout << "!!!!!!!!!!!!!!!!!!!!" <    153             G4cout << "!!!!!!!!!!!!!!!!!!!!" << G4endl;
154             G4cout << "!!! WARNING" << G4endl;    154             G4cout << "!!! WARNING" << G4endl;
155             G4cout << "G4MoleculeEncounterStep    155             G4cout << "G4MoleculeEncounterStepper::CalculateStep will return infinity "
156                 "for the reaction because the     156                 "for the reaction because the molecule "
157                 << pMoleculeA->GetName()          157                 << pMoleculeA->GetName()
158                 << " does not have any reactan    158                 << " does not have any reactants given in the reaction table."
159                 << G4endl;                        159                 << G4endl;
160             G4cout << "!!!!!!!!!!!!!!!!!!!!" <    160             G4cout << "!!!!!!!!!!!!!!!!!!!!" << G4endl;
161         }                                         161         }
162 #endif                                            162 #endif
163         return DBL_MAX;                           163         return DBL_MAX;
164     }                                             164     }
165                                                   165 
166     auto  nbReactives = (G4int)pReactantList-> << 166     G4int nbReactives = (G4int)pReactantList->size();
167                                                   167 
168     if (nbReactives == 0)                         168     if (nbReactives == 0)
169     {                                             169     {
170 #ifdef G4VERBOSE                                  170 #ifdef G4VERBOSE
171         //    DEBUG                               171         //    DEBUG
172         if (fVerbose != 0)                     << 172         if (fVerbose)
173         {                                         173         {
174             // TODO replace with the warning m    174             // TODO replace with the warning mode of G4Exception
175             G4cout << "!!!!!!!!!!!!!!!!!!!!" <    175             G4cout << "!!!!!!!!!!!!!!!!!!!!" << G4endl;
176             G4cout << "!!! WARNING" << G4endl;    176             G4cout << "!!! WARNING" << G4endl;
177             G4cout << "G4MoleculeEncounterStep    177             G4cout << "G4MoleculeEncounterStepper::CalculateStep will return infinity "
178                 "for the reaction because the     178                 "for the reaction because the molecule "
179                 << pMoleculeA->GetName()          179                 << pMoleculeA->GetName()
180                 << " does not have any reactan    180                 << " does not have any reactants given in the reaction table."
181                 << "This message can also resu    181                 << "This message can also result from a wrong implementation of the reaction table."
182                 << G4endl;                        182                 << G4endl;
183             G4cout << "!!!!!!!!!!!!!!!!!!!!" <    183             G4cout << "!!!!!!!!!!!!!!!!!!!!" << G4endl;
184         }                                         184         }
185 #endif                                            185 #endif
186         return DBL_MAX;                           186         return DBL_MAX;
187     }                                             187     }
188                                                   188 
189     fReactants = std::make_shared<vector<G4Tra << 189     fReactants.reset(new vector<G4Track*>());
190     fReactionModel->Initialise(pMolConfA, trac    190     fReactionModel->Initialise(pMolConfA, trackA);
191                                                   191 
192     //________________________________________    192     //__________________________________________________________________
193     // Start looping on possible reactants        193     // Start looping on possible reactants
194     for (G4int i = 0; i < nbReactives; i++)       194     for (G4int i = 0; i < nbReactives; i++)
195     {                                             195     {
196         auto pMoleculeB = (*pReactantList)[i];    196         auto pMoleculeB = (*pReactantList)[i];
197                                                   197 
198         //____________________________________    198         //______________________________________________________________
199         // Retrieve reaction range                199         // Retrieve reaction range
200         const G4double R = fReactionModel->Get    200         const G4double R = fReactionModel->GetReactionRadius(i);
201                                                   201 
202         //____________________________________    202         //______________________________________________________________
203         // Use KdTree algorithm to find closes    203         // Use KdTree algorithm to find closest reactants
204         G4KDTreeResultHandle resultsNearest(      204         G4KDTreeResultHandle resultsNearest(
205             G4MoleculeFinder::Instance()->Find    205             G4MoleculeFinder::Instance()->FindNearest(pMoleculeA,
206                                                   206                                                       pMoleculeB->GetMoleculeID()));
207                                                   207 
208         if (static_cast<int>(resultsNearest) = << 208         if (resultsNearest == 0) continue;
209                                                   209 
210         G4double r2 = resultsNearest->GetDista    210         G4double r2 = resultsNearest->GetDistanceSqr();
211         Utils utils(trackA, pMoleculeB);          211         Utils utils(trackA, pMoleculeB);
212                                                   212 
213         if (r2 <= R * R) // ==> Record in rang    213         if (r2 <= R * R) // ==> Record in range
214         {                                         214         {
215             // Entering in this condition may     215             // Entering in this condition may due to the fact that molecules are very close
216             // to each other                      216             // to each other
217             // Therefore, if we only take the     217             // Therefore, if we only take the nearby reactant into account, it might have already
218             // reacted. Instead, we will take     218             // reacted. Instead, we will take all possible reactants that satisfy the condition r<R
219                                                   219 
220             if (!fHasAlreadyReachedNullTime)   << 220             if (fHasAlreadyReachedNullTime == false)
221             {                                     221             {
222                 fReactants->clear();              222                 fReactants->clear();
223                 fHasAlreadyReachedNullTime = t    223                 fHasAlreadyReachedNullTime = true;
224             }                                     224             }
225                                                   225 
226             fSampledMinTimeStep = 0.;             226             fSampledMinTimeStep = 0.;
227             G4KDTreeResultHandle resultsInRang    227             G4KDTreeResultHandle resultsInRange(
228                 G4MoleculeFinder::Instance()->    228                 G4MoleculeFinder::Instance()->FindNearestInRange(pMoleculeA,
229                                                   229                                                                  pMoleculeB->GetMoleculeID(),
230                                                   230                                                                  R));
231             CheckAndRecordResults(utils,          231             CheckAndRecordResults(utils,
232 #ifdef G4VERBOSE                                  232 #ifdef G4VERBOSE
233                                   R,              233                                   R,
234 #endif                                            234 #endif
235                                   resultsInRan    235                                   resultsInRange);
236         }                                         236         }
237         else                                      237         else
238         {                                         238         {
239             G4double r = sqrt(r2);                239             G4double r = sqrt(r2);
240             G4double tempMinET = pow(r - R, 2)    240             G4double tempMinET = pow(r - R, 2) / utils.fConstant;
241             // constant = 16 * (fDA + fDB + 2*    241             // constant = 16 * (fDA + fDB + 2*sqrt(fDA*fDB))
242                                                   242 
243             if (tempMinET <= fSampledMinTimeSt    243             if (tempMinET <= fSampledMinTimeStep)
244             {                                     244             {
245                 if (fUserMinTimeStep < DBL_MAX    245                 if (fUserMinTimeStep < DBL_MAX/*IsInf(fUserMinTimeStep) == false*/
246                     && tempMinET <= fUserMinTi    246                     && tempMinET <= fUserMinTimeStep) // ==> Record in range
247                 {                                 247                 {
248                     if (fSampledMinTimeStep >     248                     if (fSampledMinTimeStep > fUserMinTimeStep)
249                     {                             249                     {
250                         fReactants->clear();      250                         fReactants->clear();
251                     }                             251                     }
252                                                   252 
253                     fSampledMinTimeStep = fUse    253                     fSampledMinTimeStep = fUserMinTimeStep;
254                                                   254 
255                     G4double range = R + sqrt(    255                     G4double range = R + sqrt(fUserMinTimeStep*utils.fConstant);
256                                                   256 
257                     G4KDTreeResultHandle resul    257                     G4KDTreeResultHandle resultsInRange(
258                         G4MoleculeFinder::Inst    258                         G4MoleculeFinder::Instance()->
259                         FindNearestInRange(pMo    259                         FindNearestInRange(pMoleculeA,
260                                            pMo    260                                            pMoleculeB->GetMoleculeID(),
261                                            ran    261                                            range));
262                                                   262 
263                     CheckAndRecordResults(util    263                     CheckAndRecordResults(utils,
264 #ifdef G4VERBOSE                                  264 #ifdef G4VERBOSE
265                                           rang    265                                           range,
266 #endif                                            266 #endif
267                                           resu    267                                           resultsInRange);
268                 }                                 268                 }
269                 else // ==> Record nearest        269                 else // ==> Record nearest
270                 {                                 270                 {
271                     if (tempMinET < fSampledMi    271                     if (tempMinET < fSampledMinTimeStep)
272                         // to avoid cases wher    272                         // to avoid cases where fSampledMinTimeStep == tempMinET
273                     {                             273                     {
274                         fSampledMinTimeStep =     274                         fSampledMinTimeStep = tempMinET;
275                         fReactants->clear();      275                         fReactants->clear();
276                     }                             276                     }
277                                                   277 
278                     CheckAndRecordResults(util    278                     CheckAndRecordResults(utils,
279 #ifdef G4VERBOSE                                  279 #ifdef G4VERBOSE
280                                           R,      280                                           R,
281 #endif                                            281 #endif
282                                           resu    282                                           resultsNearest);
283                 }                                 283                 }
284             }                                     284             }
285         }                                         285         }
286     }                                             286     }
287                                                   287 
288 #ifdef G4VERBOSE                                  288 #ifdef G4VERBOSE
289     if (fVerbose != 0)                         << 289     if (fVerbose)
290     {                                             290     {
291         G4cout << "G4MoleculeEncounterStepper:    291         G4cout << "G4MoleculeEncounterStepper::CalculateStep will finally return :"
292             << G4BestUnit(fSampledMinTimeStep,    292             << G4BestUnit(fSampledMinTimeStep, "Time") << G4endl;
293                                                   293 
294         if (fVerbose > 1)                         294         if (fVerbose > 1)
295         {                                         295         {
296             G4cout << "Selected reactants for     296             G4cout << "Selected reactants for trackA: " << pMoleculeA->GetName()
297                 << " (" << trackA.GetTrackID()    297                 << " (" << trackA.GetTrackID() << ") are: ";
298                                                   298 
299             vector<G4Track*>::iterator it;        299             vector<G4Track*>::iterator it;
300             for (it = fReactants->begin(); it     300             for (it = fReactants->begin(); it != fReactants->end(); it++)
301             {                                     301             {
302                 G4Track* trackB = *it;            302                 G4Track* trackB = *it;
303                 G4cout << GetMolecule(trackB)-    303                 G4cout << GetMolecule(trackB)->GetName() << " ("
304                     << trackB->GetTrackID() <<    304                     << trackB->GetTrackID() << ") \t ";
305             }                                     305             }
306             G4cout << G4endl;                     306             G4cout << G4endl;
307         }                                         307         }
308     }                                             308     }
309 #endif                                            309 #endif
310     return fSampledMinTimeStep;                   310     return fSampledMinTimeStep;
311 }                                                 311 }
312                                                   312 
313 void G4DNAMoleculeEncounterStepper::CheckAndRe    313 void G4DNAMoleculeEncounterStepper::CheckAndRecordResults(const Utils& utils,
314 #ifdef G4VERBOSE                                  314 #ifdef G4VERBOSE
315                                                   315                                                           const G4double R,
316 #endif                                            316 #endif
317                                                   317                                                           G4KDTreeResultHandle& results)
318 {                                                 318 {
319     if (static_cast<int>(results) == 0)        << 319     if (results == 0)
320     {                                             320     {
321 #ifdef G4VERBOSE                                  321 #ifdef G4VERBOSE
322         if (fVerbose > 1)                         322         if (fVerbose > 1)
323         {                                         323         {
324             G4cout << "No molecule " << utils.    324             G4cout << "No molecule " << utils.fpMoleculeB->GetName()
325                 << " found to react with " <<     325                 << " found to react with " << utils.fpMoleculeA->GetName()
326                 << G4endl;                        326                 << G4endl;
327         }                                         327         }
328 #endif                                            328 #endif
329         return;                                   329         return;
330     }                                             330     }
331                                                   331 
332     for (results->Rewind(); !results->End(); r    332     for (results->Rewind(); !results->End(); results->Next())
333     {                                             333     {
334         auto  reactiveB = results->GetItem<G4I << 334         G4IT* reactiveB = results->GetItem<G4IT>();
335                                                   335 
336         if (reactiveB == nullptr)              << 336         if (reactiveB == 0)
337         {                                         337         {
338             continue;                             338             continue;
339         }                                         339         }
340                                                   340 
341         G4Track *trackB = reactiveB->GetTrack(    341         G4Track *trackB = reactiveB->GetTrack();
342                                                   342 
343         if (trackB == nullptr)                 << 343         if (trackB == 0)
344         {                                         344         {
345             G4ExceptionDescription exceptionDe    345             G4ExceptionDescription exceptionDescription;
346             exceptionDescription                  346             exceptionDescription
347                 << "The reactant B found using    347                 << "The reactant B found using the MoleculeFinder does not have a valid "
348                 "track attached to it. If this    348                 "track attached to it. If this is done on purpose, please do "
349                 "not record this molecule in t    349                 "not record this molecule in the MoleculeFinder."
350                 << G4endl;                        350                 << G4endl;
351             G4Exception("G4DNAMoleculeEncounte    351             G4Exception("G4DNAMoleculeEncounterStepper::RetrieveResults",
352                         "MoleculeEncounterStep    352                         "MoleculeEncounterStepper001", FatalErrorInArgument,
353                         exceptionDescription);    353                         exceptionDescription);
354             continue;                             354             continue;
355         }                                         355         }
356                                                   356 
357         if (trackB->GetTrackStatus() != fAlive    357         if (trackB->GetTrackStatus() != fAlive)
358         {                                         358         {
359             continue;                             359             continue;
360         }                                         360         }
361                                                   361 
362         if (trackB == &utils.fpTrackA)            362         if (trackB == &utils.fpTrackA)
363         {                                         363         {
364             G4ExceptionDescription exceptionDe    364             G4ExceptionDescription exceptionDescription;
365             exceptionDescription                  365             exceptionDescription
366                 << "A track is reacting with i    366                 << "A track is reacting with itself (which is impossible) ie fpTrackA == trackB"
367                 << G4endl;                        367                 << G4endl;
368             exceptionDescription << "Molecule     368             exceptionDescription << "Molecule A (and B) is of type : "
369                 << utils.fpMoleculeA->GetName(    369                 << utils.fpMoleculeA->GetName() << " with trackID : "
370                 << utils.fpTrackA.GetTrackID()    370                 << utils.fpTrackA.GetTrackID() << G4endl;
371                                                   371 
372             G4Exception("G4DNAMoleculeEncounte    372             G4Exception("G4DNAMoleculeEncounterStepper::RetrieveResults",
373                         "MoleculeEncounterStep    373                         "MoleculeEncounterStepper003", FatalErrorInArgument,
374                         exceptionDescription);    374                         exceptionDescription);
375                                                   375 
376         }                                         376         }
377                                                   377 
378         if (fabs(trackB->GetGlobalTime() - uti    378         if (fabs(trackB->GetGlobalTime() - utils.fpTrackA.GetGlobalTime())
379         > utils.fpTrackA.GetGlobalTime() * (1     379         > utils.fpTrackA.GetGlobalTime() * (1 - 1 / 100))
380         {                                         380         {
381             // DEBUG                              381             // DEBUG
382             G4ExceptionDescription exceptionDe    382             G4ExceptionDescription exceptionDescription;
383             exceptionDescription                  383             exceptionDescription
384                 << "The interacting tracks are    384                 << "The interacting tracks are not synchronized in time" << G4endl;
385             exceptionDescription                  385             exceptionDescription
386                 << "trackB->GetGlobalTime() !=    386                 << "trackB->GetGlobalTime() != fpTrackA.GetGlobalTime()" << G4endl;
387                                                   387 
388             exceptionDescription << "fpTrackA     388             exceptionDescription << "fpTrackA : trackID : " << utils.fpTrackA.GetTrackID()
389                 << "\t Name :" << utils.fpMole    389                 << "\t Name :" << utils.fpMoleculeA->GetName()
390                 << "\t fpTrackA->GetGlobalTime    390                 << "\t fpTrackA->GetGlobalTime() = "
391                 << G4BestUnit(utils.fpTrackA.G    391                 << G4BestUnit(utils.fpTrackA.GetGlobalTime(), "Time") << G4endl;
392                                                   392 
393             exceptionDescription << "trackB :     393             exceptionDescription << "trackB : trackID : " << trackB->GetTrackID()
394                 << "\t Name :" << utils.fpMole    394                 << "\t Name :" << utils.fpMoleculeB->GetName()
395                 << "\t trackB->GetGlobalTime()    395                 << "\t trackB->GetGlobalTime() = "
396                 << G4BestUnit(trackB->GetGloba    396                 << G4BestUnit(trackB->GetGlobalTime(), "Time") << G4endl;
397                                                   397 
398             G4Exception("G4DNAMoleculeEncounte    398             G4Exception("G4DNAMoleculeEncounterStepper::RetrieveResults",
399                         "MoleculeEncounterStep    399                         "MoleculeEncounterStepper004", FatalErrorInArgument,
400                         exceptionDescription);    400                         exceptionDescription);
401         }                                         401         }
402                                                   402 
403 #ifdef G4VERBOSE                                  403 #ifdef G4VERBOSE
404         if (fVerbose > 1)                         404         if (fVerbose > 1)
405         {                                         405         {
406                                                   406 
407             G4double r2 = results->GetDistance    407             G4double r2 = results->GetDistanceSqr();
408             G4cout << "\t ********************    408             G4cout << "\t ************************************************** " << G4endl;
409             G4cout << "\t Reaction between "      409             G4cout << "\t Reaction between "
410                 << utils.fpMoleculeA->GetName(    410                 << utils.fpMoleculeA->GetName() << " (" << utils.fpTrackA.GetTrackID() << ") "
411                 << " & " << utils.fpMoleculeB-    411                 << " & " << utils.fpMoleculeB->GetName() << " (" << trackB->GetTrackID() << "), "
412                 << "Interaction Range = "         412                 << "Interaction Range = "
413                 << G4BestUnit(R, "Length") <<     413                 << G4BestUnit(R, "Length") << G4endl;
414             G4cout << "\t Real distance betwee    414             G4cout << "\t Real distance between reactants  = "
415                 << G4BestUnit((utils.fpTrackA.    415                 << G4BestUnit((utils.fpTrackA.GetPosition() - trackB->GetPosition()).mag(), "Length") << G4endl;
416             G4cout << "\t Distance between rea    416             G4cout << "\t Distance between reactants calculated by nearest neighbor algorithm = "
417                 << G4BestUnit(sqrt(r2), "Lengt    417                 << G4BestUnit(sqrt(r2), "Length") << G4endl;
418                                                   418 
419         }                                         419         }
420 #endif                                            420 #endif
421                                                   421 
422         fReactants->push_back(trackB);            422         fReactants->push_back(trackB);
423     }                                             423     }
424 }                                                 424 }
425                                                   425 
426 void G4DNAMoleculeEncounterStepper::SetReactio    426 void G4DNAMoleculeEncounterStepper::SetReactionModel(G4VDNAReactionModel* pReactionModel)
427 {                                                 427 {
428     fReactionModel = pReactionModel;              428     fReactionModel = pReactionModel;
429 }                                                 429 }
430                                                   430 
431 G4VDNAReactionModel* G4DNAMoleculeEncounterSte    431 G4VDNAReactionModel* G4DNAMoleculeEncounterStepper::GetReactionModel()
432 {                                                 432 {
433     return fReactionModel;                        433     return fReactionModel;
434 }                                                 434 }
435                                                   435 
436 void G4DNAMoleculeEncounterStepper::SetVerbose    436 void G4DNAMoleculeEncounterStepper::SetVerbose(int flag)
437 {                                                 437 {
438     fVerbose = flag;                              438     fVerbose = flag;
439 }                                                 439 }
440                                                   440 
441 G4double G4DNAMoleculeEncounterStepper::Calcul    441 G4double G4DNAMoleculeEncounterStepper::CalculateMinTimeStep(G4double /*currentGlobalTime*/, G4double definedMinTimeStep){
442                                                   442 
443     G4double fTSTimeStep = DBL_MAX;               443     G4double fTSTimeStep = DBL_MAX;
444                                                   444 
445     for (auto pTrack : *fpTrackContainer->GetM    445     for (auto pTrack : *fpTrackContainer->GetMainList())
446     {                                             446     {
447         if (pTrack == nullptr)                    447         if (pTrack == nullptr)
448         {                                         448         {
449             G4ExceptionDescription exceptionDe    449             G4ExceptionDescription exceptionDescription;
450             exceptionDescription << "No track     450             exceptionDescription << "No track found.";
451             G4Exception("G4Scheduler::Calculat    451             G4Exception("G4Scheduler::CalculateMinStep", "ITScheduler006",
452                         FatalErrorInArgument,     452                         FatalErrorInArgument, exceptionDescription);
453             continue;                             453             continue;
454         }                                         454         }
455                                                   455 
456         G4TrackStatus trackStatus = pTrack->Ge    456         G4TrackStatus trackStatus = pTrack->GetTrackStatus();
457         if (trackStatus == fStopAndKill || tra    457         if (trackStatus == fStopAndKill || trackStatus == fStopButAlive)
458         {                                         458         {
459             continue;                             459             continue;
460         }                                         460         }
461                                                   461 
462         G4double sampledMinTimeStep = Calculat    462         G4double sampledMinTimeStep = CalculateStep(*pTrack, definedMinTimeStep);
463         G4TrackVectorHandle reactants = GetRea    463         G4TrackVectorHandle reactants = GetReactants();
464                                                   464 
465         if (sampledMinTimeStep < fTSTimeStep)     465         if (sampledMinTimeStep < fTSTimeStep)
466         {                                         466         {
467             fTSTimeStep = sampledMinTimeStep;     467             fTSTimeStep = sampledMinTimeStep;
468             fReactionSet->CleanAllReaction();     468             fReactionSet->CleanAllReaction();
469             if (reactants)                        469             if (reactants)
470             {                                     470             {
471                 fReactionSet->AddReactions(fTS    471                 fReactionSet->AddReactions(fTSTimeStep,
472                                            con    472                                            const_cast<G4Track*>(pTrack),
473                                            std << 473                                            reactants);
474                 ResetReactants();                 474                 ResetReactants();
475             }                                     475             }
476          }                                        476          }
477          else if (fTSTimeStep == sampledMinTim    477          else if (fTSTimeStep == sampledMinTimeStep && bool(reactants))
478          {                                        478          {
479              fReactionSet->AddReactions(fTSTim    479              fReactionSet->AddReactions(fTSTimeStep,
480                                         const_    480                                         const_cast<G4Track*>(pTrack),
481                                         std::m << 481                                         reactants);
482              ResetReactants();                    482              ResetReactants();
483          }                                        483          }
484          else if (reactants)                      484          else if (reactants)
485          {                                        485          {
486              ResetReactants();                    486              ResetReactants();
487          }                                        487          }
488     }                                             488     }
489                                                   489 
490     return fTSTimeStep;                           490     return fTSTimeStep;
491 }                                                 491 }
492                                                   492