| Geant4 Cross Reference |
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
40 #include "G4DNAMolecularReactionTable.hh" 40 #include "G4DNAMolecularReactionTable.hh"
41 #include "G4H2O.hh" 41 #include "G4H2O.hh"
42 #include "G4MolecularConfiguration.hh" 42 #include "G4MolecularConfiguration.hh"
43 #include "G4MoleculeFinder.hh" 43 #include "G4MoleculeFinder.hh"
44 #include "G4UnitsTable.hh" 44 #include "G4UnitsTable.hh"
45 #include "G4VDNAReactionModel.hh" 45 #include "G4VDNAReactionModel.hh"
46 #include "G4memory.hh" 46 #include "G4memory.hh"
47 47
48 #include <memory> 48 #include <memory>
49 49
50 using namespace std; 50 using namespace std;
51 using namespace CLHEP; 51 using namespace CLHEP;
52 52
53 //#define DEBUG_MEM 53 //#define DEBUG_MEM
54 54
55 #ifdef DEBUG_MEM 55 #ifdef DEBUG_MEM
56 #include "G4MemStat.hh" 56 #include "G4MemStat.hh"
57 using namespace G4MemStat; 57 using namespace G4MemStat;
58 #endif 58 #endif
59 59
60 G4DNAMoleculeEncounterStepper::Utils::Utils(co 60 G4DNAMoleculeEncounterStepper::Utils::Utils(const G4Track& tA,
61 co 61 const G4MolecularConfiguration* pMoleculeB)
62 : fpTrackA(tA) 62 : fpTrackA(tA)
63 , fpMoleculeB(pMoleculeB) 63 , fpMoleculeB(pMoleculeB)
64 { 64 {
65 fpMoleculeA = GetMolecule(tA); 65 fpMoleculeA = GetMolecule(tA);
66 fDA = fpMoleculeA->GetDiffusionCoefficient 66 fDA = fpMoleculeA->GetDiffusionCoefficient();
67 fDB = fpMoleculeB->GetDiffusionCoefficient 67 fDB = fpMoleculeB->GetDiffusionCoefficient();
68 fConstant = 8 * (fDA + fDB + 2 * sqrt(fDA 68 fConstant = 8 * (fDA + fDB + 2 * sqrt(fDA * fDB));
69 } 69 }
70 70
71 G4DNAMoleculeEncounterStepper::G4DNAMoleculeEn 71 G4DNAMoleculeEncounterStepper::G4DNAMoleculeEncounterStepper()
72 : 72 :
73 fMolecularReactionTable(reference_cast<co 73 fMolecularReactionTable(reference_cast<const G4DNAMolecularReactionTable*>(fpReactionTable))
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 != 0)
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 == nullptr)
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 auto 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 != 0)
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 = std::make_shared<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 (static_cast<int>(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)
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 != 0)
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 (static_cast<int>(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 auto reactiveB = results->GetItem<G4IT>();
335 335
336 if (reactiveB == nullptr) 336 if (reactiveB == nullptr)
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 == nullptr)
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