| 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 /* 27 /*
28 * G4DNAIRTMoleculeEncounterStepper.cc 28 * G4DNAIRTMoleculeEncounterStepper.cc
29 * 29 *
30 * Created on: Jul 23, 2019 30 * Created on: Jul 23, 2019
31 * Author: W. G. Shin 31 * Author: W. G. Shin
32 * J. Ramos-Mendez and B. Faddego 32 * J. Ramos-Mendez and B. Faddegon
33 */ 33 */
34 34
35 #include "G4DNAIRTMoleculeEncounterStepper.hh" 35 #include "G4DNAIRTMoleculeEncounterStepper.hh"
36 36
37 #include "G4DNAMolecularReactionTable.hh" 37 #include "G4DNAMolecularReactionTable.hh"
38 #include "G4H2O.hh" 38 #include "G4H2O.hh"
39 #include "G4ITReaction.hh" 39 #include "G4ITReaction.hh"
40 #include "G4MolecularConfiguration.hh" 40 #include "G4MolecularConfiguration.hh"
41 #include "G4MoleculeFinder.hh" 41 #include "G4MoleculeFinder.hh"
42 #include "G4Scheduler.hh" 42 #include "G4Scheduler.hh"
43 #include "G4UnitsTable.hh" 43 #include "G4UnitsTable.hh"
44 #include "G4VDNAReactionModel.hh" 44 #include "G4VDNAReactionModel.hh"
45 #include "G4memory.hh" 45 #include "G4memory.hh"
46 46
47 #include <memory> 47 #include <memory>
48 48
49 using namespace std; 49 using namespace std;
50 using namespace CLHEP; 50 using namespace CLHEP;
51 51
52 //#define DEBUG_MEM 52 //#define DEBUG_MEM
53 53
54 #ifdef DEBUG_MEM 54 #ifdef DEBUG_MEM
55 #include "G4MemStat.hh" 55 #include "G4MemStat.hh"
56 using namespace G4MemStat; 56 using namespace G4MemStat;
57 #endif 57 #endif
58 58
59 G4DNAIRTMoleculeEncounterStepper::Utils::Utils 59 G4DNAIRTMoleculeEncounterStepper::Utils::Utils(const G4Track& tA,
60 co 60 const G4MolecularConfiguration* pMoleculeB)
61 : fpTrackA(tA) 61 : fpTrackA(tA)
62 , fpMoleculeB(pMoleculeB) 62 , fpMoleculeB(pMoleculeB)
63 { 63 {
64 fpMoleculeA = GetMolecule(tA); 64 fpMoleculeA = GetMolecule(tA);
65 fDA = fpMoleculeA->GetDiffusionCoefficient 65 fDA = fpMoleculeA->GetDiffusionCoefficient();
66 fDB = fpMoleculeB->GetDiffusionCoefficient 66 fDB = fpMoleculeB->GetDiffusionCoefficient();
67 fConstant = 8 * (fDA + fDB + 2 * sqrt(fDA 67 fConstant = 8 * (fDA + fDB + 2 * sqrt(fDA * fDB));
68 } 68 }
69 69
70 G4DNAIRTMoleculeEncounterStepper::G4DNAIRTMole 70 G4DNAIRTMoleculeEncounterStepper::G4DNAIRTMoleculeEncounterStepper()
71 : 71 :
72 fMolecularReactionTable(reference_cast<co 72 fMolecularReactionTable(reference_cast<const G4DNAMolecularReactionTable*>(fpReactionTable))
73 { 73 {
74 fpTrackContainer = G4ITTrackHolder::Instance 74 fpTrackContainer = G4ITTrackHolder::Instance();
75 fReactionSet = G4ITReactionSet::Instance(); 75 fReactionSet = G4ITReactionSet::Instance();
76 } 76 }
77 77
78 G4DNAIRTMoleculeEncounterStepper::~G4DNAIRTMol 78 G4DNAIRTMoleculeEncounterStepper::~G4DNAIRTMoleculeEncounterStepper() = default;
79 79
80 void G4DNAIRTMoleculeEncounterStepper::Prepare 80 void G4DNAIRTMoleculeEncounterStepper::Prepare()
81 { 81 {
82 fSampledMinTimeStep = DBL_MAX; 82 fSampledMinTimeStep = DBL_MAX;
83 if(G4Scheduler::Instance()->GetGlobalTime( 83 if(G4Scheduler::Instance()->GetGlobalTime() == G4Scheduler::Instance()->GetStartTime()){
84 G4VITTimeStepComputer::Prepare(); 84 G4VITTimeStepComputer::Prepare();
85 G4MoleculeFinder::Instance()->UpdatePo 85 G4MoleculeFinder::Instance()->UpdatePositionMap();
86 } 86 }
87 } 87 }
88 88
89 void G4DNAIRTMoleculeEncounterStepper::Initial 89 void G4DNAIRTMoleculeEncounterStepper::InitializeForNewTrack()
90 { 90 {
91 if (fReactants) 91 if (fReactants)
92 { 92 {
93 fReactants.reset(); 93 fReactants.reset();
94 } 94 }
95 fSampledMinTimeStep = DBL_MAX; 95 fSampledMinTimeStep = DBL_MAX;
96 fHasAlreadyReachedNullTime = false; 96 fHasAlreadyReachedNullTime = false;
97 } 97 }
98 98
99 template<typename T> 99 template<typename T>
100 inline bool IsInf(T value) 100 inline bool IsInf(T value)
101 { 101 {
102 return std::numeric_limits<T>::has_infinit 102 return std::numeric_limits<T>::has_infinity
103 && value == std::numeric_limits<T>::in 103 && value == std::numeric_limits<T>::infinity();
104 } 104 }
105 105
106 G4double 106 G4double
107 G4DNAIRTMoleculeEncounterStepper::CalculateSte 107 G4DNAIRTMoleculeEncounterStepper::CalculateStep(const G4Track& trackA,
108 c 108 const G4double& userMinTimeStep)
109 { 109 {
110 110
111 auto pMoleculeA = GetMolecule(trackA); 111 auto pMoleculeA = GetMolecule(trackA);
112 InitializeForNewTrack(); 112 InitializeForNewTrack();
113 fUserMinTimeStep = userMinTimeStep; 113 fUserMinTimeStep = userMinTimeStep;
114 114
115 #ifdef G4VERBOSE 115 #ifdef G4VERBOSE
116 if (fVerbose != 0) 116 if (fVerbose != 0)
117 { 117 {
118 G4cout 118 G4cout
119 << "______________________________ 119 << "_______________________________________________________________________"
120 << G4endl; 120 << G4endl;
121 G4cout << "G4DNAMoleculeEncounterStepp 121 G4cout << "G4DNAMoleculeEncounterStepper::CalculateStep" << G4endl;
122 G4cout << "Check done for molecule : " 122 G4cout << "Check done for molecule : " << pMoleculeA->GetName()
123 << " (" << trackA.GetTrackID() << 123 << " (" << trackA.GetTrackID() << ") "
124 << G4endl; 124 << G4endl;
125 } 125 }
126 #endif 126 #endif
127 127
128 //________________________________________ 128 //__________________________________________________________________
129 // Retrieve general informations for makin 129 // Retrieve general informations for making reactions
130 auto pMolConfA = pMoleculeA->GetMolecularC 130 auto pMolConfA = pMoleculeA->GetMolecularConfiguration();
131 131
132 const auto pReactantList = fMolecularReact 132 const auto pReactantList = fMolecularReactionTable->CanReactWith(pMolConfA);
133 133
134 if (pReactantList == nullptr) 134 if (pReactantList == nullptr)
135 { 135 {
136 #ifdef G4VERBOSE 136 #ifdef G4VERBOSE
137 // DEBUG 137 // DEBUG
138 if (fVerbose > 1) 138 if (fVerbose > 1)
139 { 139 {
140 G4cout << "!!!!!!!!!!!!!!!!!!!!" < 140 G4cout << "!!!!!!!!!!!!!!!!!!!!" << G4endl;
141 G4cout << "!!! WARNING" << G4endl; 141 G4cout << "!!! WARNING" << G4endl;
142 G4cout << "G4MoleculeEncounterStep 142 G4cout << "G4MoleculeEncounterStepper::CalculateStep will return infinity "
143 "for the reaction because the 143 "for the reaction because the molecule "
144 << pMoleculeA->GetName() 144 << pMoleculeA->GetName()
145 << " does not have any reactan 145 << " does not have any reactants given in the reaction table."
146 << G4endl; 146 << G4endl;
147 G4cout << "!!!!!!!!!!!!!!!!!!!!" < 147 G4cout << "!!!!!!!!!!!!!!!!!!!!" << G4endl;
148 } 148 }
149 #endif 149 #endif
150 return DBL_MAX; 150 return DBL_MAX;
151 } 151 }
152 152
153 auto nbReactives = (G4int)pReactantList-> 153 auto nbReactives = (G4int)pReactantList->size();
154 154
155 if (nbReactives == 0) 155 if (nbReactives == 0)
156 { 156 {
157 #ifdef G4VERBOSE 157 #ifdef G4VERBOSE
158 // DEBUG 158 // DEBUG
159 if (fVerbose != 0) 159 if (fVerbose != 0)
160 { 160 {
161 // TODO replace with the warning m 161 // TODO replace with the warning mode of G4Exception
162 G4cout << "!!!!!!!!!!!!!!!!!!!!" < 162 G4cout << "!!!!!!!!!!!!!!!!!!!!" << G4endl;
163 G4cout << "!!! WARNING" << G4endl; 163 G4cout << "!!! WARNING" << G4endl;
164 G4cout << "G4MoleculeEncounterStep 164 G4cout << "G4MoleculeEncounterStepper::CalculateStep will return infinity "
165 "for the reaction because the 165 "for the reaction because the molecule "
166 << pMoleculeA->GetName() 166 << pMoleculeA->GetName()
167 << " does not have any reactan 167 << " does not have any reactants given in the reaction table."
168 << "This message can also resu 168 << "This message can also result from a wrong implementation of the reaction table."
169 << G4endl; 169 << G4endl;
170 G4cout << "!!!!!!!!!!!!!!!!!!!!" < 170 G4cout << "!!!!!!!!!!!!!!!!!!!!" << G4endl;
171 } 171 }
172 #endif 172 #endif
173 return DBL_MAX; 173 return DBL_MAX;
174 } 174 }
175 175
176 fReactants = std::make_shared<vector<G4Tra 176 fReactants = std::make_shared<vector<G4Track*>>();
177 fReactionModel->Initialise(pMolConfA, trac 177 fReactionModel->Initialise(pMolConfA, trackA);
178 178
179 //________________________________________ 179 //__________________________________________________________________
180 // Start looping on possible reactants 180 // Start looping on possible reactants
181 for (G4int i = 0; i < nbReactives; ++i) 181 for (G4int i = 0; i < nbReactives; ++i)
182 { 182 {
183 auto pMoleculeB = (*pReactantList)[i]; 183 auto pMoleculeB = (*pReactantList)[i];
184 184
185 //____________________________________ 185 //______________________________________________________________
186 // Retrieve reaction range 186 // Retrieve reaction range
187 const G4double R = fReactionModel->Get 187 const G4double R = fReactionModel->GetReactionRadius(i);
188 188
189 //____________________________________ 189 //______________________________________________________________
190 // Use KdTree algorithm to find closes 190 // Use KdTree algorithm to find closest reactants
191 G4KDTreeResultHandle resultsNearest( 191 G4KDTreeResultHandle resultsNearest(
192 G4MoleculeFinder::Instance()->Find 192 G4MoleculeFinder::Instance()->FindNearest(pMoleculeA,
193 193 pMoleculeB->GetMoleculeID()));
194 194
195 if (static_cast<int>(resultsNearest) = 195 if (static_cast<int>(resultsNearest) == 0) continue;
196 196
197 G4double r2 = resultsNearest->GetDista 197 G4double r2 = resultsNearest->GetDistanceSqr();
198 Utils utils(trackA, pMoleculeB); 198 Utils utils(trackA, pMoleculeB);
199 199
200 if (r2 <= R * R) // ==> Record in rang 200 if (r2 <= R * R) // ==> Record in range
201 { 201 {
202 // Entering in this condition may 202 // Entering in this condition may due to the fact that molecules are very close
203 // to each other 203 // to each other
204 // Therefore, if we only take the 204 // Therefore, if we only take the nearby reactant into account, it might have already
205 // reacted. Instead, we will take 205 // reacted. Instead, we will take all possible reactants that satisfy the condition r<R
206 206
207 if (!fHasAlreadyReachedNullTime) 207 if (!fHasAlreadyReachedNullTime)
208 { 208 {
209 fReactants->clear(); 209 fReactants->clear();
210 fHasAlreadyReachedNullTime = t 210 fHasAlreadyReachedNullTime = true;
211 } 211 }
212 212
213 fSampledMinTimeStep = 0.; 213 fSampledMinTimeStep = 0.;
214 G4KDTreeResultHandle resultsInRang 214 G4KDTreeResultHandle resultsInRange(
215 G4MoleculeFinder::Instance()-> 215 G4MoleculeFinder::Instance()->FindNearestInRange(pMoleculeA,
216 216 pMoleculeB->GetMoleculeID(),
217 217 R));
218 CheckAndRecordResults(utils, 218 CheckAndRecordResults(utils,
219 #ifdef G4VERBOSE 219 #ifdef G4VERBOSE
220 R, 220 R,
221 #endif 221 #endif
222 resultsInRan 222 resultsInRange);
223 } 223 }
224 else 224 else
225 { 225 {
226 G4double r = sqrt(r2); 226 G4double r = sqrt(r2);
227 G4double tempMinET = pow(r - R, 2) 227 G4double tempMinET = pow(r - R, 2) / utils.fConstant;
228 // constant = 16 * (fDA + fDB + 2* 228 // constant = 16 * (fDA + fDB + 2*sqrt(fDA*fDB))
229 229
230 if (tempMinET <= fSampledMinTimeSt 230 if (tempMinET <= fSampledMinTimeStep)
231 { 231 {
232 if (fUserMinTimeStep < DBL_MAX 232 if (fUserMinTimeStep < DBL_MAX/*IsInf(fUserMinTimeStep) == false*/
233 && tempMinET <= fUserMinTi 233 && tempMinET <= fUserMinTimeStep) // ==> Record in range
234 { 234 {
235 if (fSampledMinTimeStep > 235 if (fSampledMinTimeStep > fUserMinTimeStep)
236 { 236 {
237 fReactants->clear(); 237 fReactants->clear();
238 } 238 }
239 239
240 fSampledMinTimeStep = fUse 240 fSampledMinTimeStep = fUserMinTimeStep;
241 241
242 G4double range = R + sqrt( 242 G4double range = R + sqrt(fUserMinTimeStep*utils.fConstant);
243 243
244 G4KDTreeResultHandle resul 244 G4KDTreeResultHandle resultsInRange(
245 G4MoleculeFinder::Inst 245 G4MoleculeFinder::Instance()->
246 FindNearestInRange(pMo 246 FindNearestInRange(pMoleculeA,
247 pMo 247 pMoleculeB->GetMoleculeID(),
248 ran 248 range));
249 249
250 CheckAndRecordResults(util 250 CheckAndRecordResults(utils,
251 #ifdef G4VERBOSE 251 #ifdef G4VERBOSE
252 rang 252 range,
253 #endif 253 #endif
254 resu 254 resultsInRange);
255 } 255 }
256 else // ==> Record nearest 256 else // ==> Record nearest
257 { 257 {
258 if (tempMinET < fSampledMi 258 if (tempMinET < fSampledMinTimeStep)
259 // to avoid cases wher 259 // to avoid cases where fSampledMinTimeStep == tempMinET
260 { 260 {
261 fSampledMinTimeStep = 261 fSampledMinTimeStep = tempMinET;
262 fReactants->clear(); 262 fReactants->clear();
263 } 263 }
264 264
265 CheckAndRecordResults(util 265 CheckAndRecordResults(utils,
266 #ifdef G4VERBOSE 266 #ifdef G4VERBOSE
267 R, 267 R,
268 #endif 268 #endif
269 resu 269 resultsNearest);
270 } 270 }
271 } 271 }
272 } 272 }
273 } 273 }
274 274
275 #ifdef G4VERBOSE 275 #ifdef G4VERBOSE
276 if (fVerbose != 0) 276 if (fVerbose != 0)
277 { 277 {
278 G4cout << "G4MoleculeEncounterStepper: 278 G4cout << "G4MoleculeEncounterStepper::CalculateStep will finally return :"
279 << G4BestUnit(fSampledMinTimeStep, 279 << G4BestUnit(fSampledMinTimeStep, "Time") << G4endl;
280 280
281 if (fVerbose > 1) 281 if (fVerbose > 1)
282 { 282 {
283 G4cout << "Selected reactants for 283 G4cout << "Selected reactants for trackA: " << pMoleculeA->GetName()
284 << " (" << trackA.GetTrackID() 284 << " (" << trackA.GetTrackID() << ") are: ";
285 285
286 vector<G4Track*>::iterator it; 286 vector<G4Track*>::iterator it;
287 for (it = fReactants->begin(); it 287 for (it = fReactants->begin(); it != fReactants->end(); it++)
288 { 288 {
289 G4Track* trackB = *it; 289 G4Track* trackB = *it;
290 G4cout << GetMolecule(trackB)- 290 G4cout << GetMolecule(trackB)->GetName() << " ("
291 << trackB->GetTrackID() << 291 << trackB->GetTrackID() << ") \t ";
292 } 292 }
293 G4cout << G4endl; 293 G4cout << G4endl;
294 } 294 }
295 } 295 }
296 #endif 296 #endif
297 return fSampledMinTimeStep; 297 return fSampledMinTimeStep;
298 } 298 }
299 299
300 300
301 301
302 302
303 void G4DNAIRTMoleculeEncounterStepper::CheckAn 303 void G4DNAIRTMoleculeEncounterStepper::CheckAndRecordResults(const Utils& utils,
304 #ifdef G4VERBOSE 304 #ifdef G4VERBOSE
305 305 const G4double R,
306 #endif 306 #endif
307 307 G4KDTreeResultHandle& results)
308 { 308 {
309 if (static_cast<int>(results) == 0) 309 if (static_cast<int>(results) == 0)
310 { 310 {
311 #ifdef G4VERBOSE 311 #ifdef G4VERBOSE
312 if (fVerbose > 1) 312 if (fVerbose > 1)
313 { 313 {
314 G4cout << "No molecule " << utils. 314 G4cout << "No molecule " << utils.fpMoleculeB->GetName()
315 << " found to react with " << 315 << " found to react with " << utils.fpMoleculeA->GetName()
316 << G4endl; 316 << G4endl;
317 } 317 }
318 #endif 318 #endif
319 return; 319 return;
320 } 320 }
321 321
322 for (results->Rewind(); !results->End(); r 322 for (results->Rewind(); !results->End(); results->Next())
323 { 323 {
324 auto reactiveB = results->GetItem<G4I 324 auto reactiveB = results->GetItem<G4IT>();
325 325
326 if (reactiveB == nullptr) 326 if (reactiveB == nullptr)
327 { 327 {
328 continue; 328 continue;
329 } 329 }
330 330
331 G4Track *trackB = reactiveB->GetTrack( 331 G4Track *trackB = reactiveB->GetTrack();
332 332
333 if (trackB == nullptr) 333 if (trackB == nullptr)
334 { 334 {
335 G4ExceptionDescription exceptionDe 335 G4ExceptionDescription exceptionDescription;
336 exceptionDescription 336 exceptionDescription
337 << "The reactant B found using 337 << "The reactant B found using the MoleculeFinder does not have a valid "
338 "track attached to it. If this 338 "track attached to it. If this is done on purpose, please do "
339 "not record this molecule in t 339 "not record this molecule in the MoleculeFinder."
340 << G4endl; 340 << G4endl;
341 G4Exception("G4DNAMoleculeEncounte 341 G4Exception("G4DNAMoleculeEncounterStepper::RetrieveResults",
342 "MoleculeEncounterStep 342 "MoleculeEncounterStepper001", FatalErrorInArgument,
343 exceptionDescription); 343 exceptionDescription);
344 continue; 344 continue;
345 } 345 }
346 346
347 if (trackB->GetTrackStatus() != fAlive 347 if (trackB->GetTrackStatus() != fAlive)
348 { 348 {
349 continue; 349 continue;
350 } 350 }
351 351
352 if (trackB == &utils.fpTrackA) 352 if (trackB == &utils.fpTrackA)
353 { 353 {
354 G4ExceptionDescription exceptionDe 354 G4ExceptionDescription exceptionDescription;
355 exceptionDescription 355 exceptionDescription
356 << "A track is reacting with i 356 << "A track is reacting with itself (which is impossible) ie fpTrackA == trackB"
357 << G4endl; 357 << G4endl;
358 exceptionDescription << "Molecule 358 exceptionDescription << "Molecule A (and B) is of type : "
359 << utils.fpMoleculeA->GetName( 359 << utils.fpMoleculeA->GetName() << " with trackID : "
360 << utils.fpTrackA.GetTrackID() 360 << utils.fpTrackA.GetTrackID() << G4endl;
361 361
362 G4Exception("G4DNAMoleculeEncounte 362 G4Exception("G4DNAMoleculeEncounterStepper::RetrieveResults",
363 "MoleculeEncounterStep 363 "MoleculeEncounterStepper003", FatalErrorInArgument,
364 exceptionDescription); 364 exceptionDescription);
365 365
366 } 366 }
367 367
368 if (fabs(trackB->GetGlobalTime() - uti 368 if (fabs(trackB->GetGlobalTime() - utils.fpTrackA.GetGlobalTime())
369 > utils.fpTrackA.GetGlobalTime() * (1 369 > utils.fpTrackA.GetGlobalTime() * (1 - 1 / 100))
370 { 370 {
371 // DEBUG 371 // DEBUG
372 G4ExceptionDescription exceptionDe 372 G4ExceptionDescription exceptionDescription;
373 exceptionDescription 373 exceptionDescription
374 << "The interacting tracks are 374 << "The interacting tracks are not synchronized in time" << G4endl;
375 exceptionDescription 375 exceptionDescription
376 << "trackB->GetGlobalTime() != 376 << "trackB->GetGlobalTime() != fpTrackA.GetGlobalTime()" << G4endl;
377 377
378 exceptionDescription << "fpTrackA 378 exceptionDescription << "fpTrackA : trackID : " << utils.fpTrackA.GetTrackID()
379 << "\t Name :" << utils.fpMole 379 << "\t Name :" << utils.fpMoleculeA->GetName()
380 << "\t fpTrackA->GetGlobalTime 380 << "\t fpTrackA->GetGlobalTime() = "
381 << G4BestUnit(utils.fpTrackA.G 381 << G4BestUnit(utils.fpTrackA.GetGlobalTime(), "Time") << G4endl;
382 382
383 exceptionDescription << "trackB : 383 exceptionDescription << "trackB : trackID : " << trackB->GetTrackID()
384 << "\t Name :" << utils.fpMole 384 << "\t Name :" << utils.fpMoleculeB->GetName()
385 << "\t trackB->GetGlobalTime() 385 << "\t trackB->GetGlobalTime() = "
386 << G4BestUnit(trackB->GetGloba 386 << G4BestUnit(trackB->GetGlobalTime(), "Time") << G4endl;
387 387
388 G4Exception("G4DNAMoleculeEncounte 388 G4Exception("G4DNAMoleculeEncounterStepper::RetrieveResults",
389 "MoleculeEncounterStep 389 "MoleculeEncounterStepper004", FatalErrorInArgument,
390 exceptionDescription); 390 exceptionDescription);
391 } 391 }
392 392
393 #ifdef G4VERBOSE 393 #ifdef G4VERBOSE
394 if (fVerbose > 1) 394 if (fVerbose > 1)
395 { 395 {
396 396
397 G4double r2 = results->GetDistance 397 G4double r2 = results->GetDistanceSqr();
398 G4cout << "\t ******************** 398 G4cout << "\t ************************************************** " << G4endl;
399 G4cout << "\t Reaction between " 399 G4cout << "\t Reaction between "
400 << utils.fpMoleculeA->GetName( 400 << utils.fpMoleculeA->GetName() << " (" << utils.fpTrackA.GetTrackID() << ") "
401 << " & " << utils.fpMoleculeB- 401 << " & " << utils.fpMoleculeB->GetName() << " (" << trackB->GetTrackID() << "), "
402 << "Interaction Range = " 402 << "Interaction Range = "
403 << G4BestUnit(R, "Length") << 403 << G4BestUnit(R, "Length") << G4endl;
404 G4cout << "\t Real distance betwee 404 G4cout << "\t Real distance between reactants = "
405 << G4BestUnit((utils.fpTrackA. 405 << G4BestUnit((utils.fpTrackA.GetPosition() - trackB->GetPosition()).mag(), "Length") << G4endl;
406 G4cout << "\t Distance between rea 406 G4cout << "\t Distance between reactants calculated by nearest neighbor algorithm = "
407 << G4BestUnit(sqrt(r2), "Lengt 407 << G4BestUnit(sqrt(r2), "Length") << G4endl;
408 408
409 } 409 }
410 #endif 410 #endif
411 411
412 fReactants->push_back(trackB); 412 fReactants->push_back(trackB);
413 } 413 }
414 } 414 }
415 415
416 void G4DNAIRTMoleculeEncounterStepper::SetReac 416 void G4DNAIRTMoleculeEncounterStepper::SetReactionModel(G4VDNAReactionModel* pReactionModel)
417 { 417 {
418 fReactionModel = pReactionModel; 418 fReactionModel = pReactionModel;
419 } 419 }
420 420
421 G4VDNAReactionModel* G4DNAIRTMoleculeEncounter 421 G4VDNAReactionModel* G4DNAIRTMoleculeEncounterStepper::GetReactionModel()
422 { 422 {
423 return fReactionModel; 423 return fReactionModel;
424 } 424 }
425 425
426 void G4DNAIRTMoleculeEncounterStepper::SetVerb 426 void G4DNAIRTMoleculeEncounterStepper::SetVerbose(int flag)
427 { 427 {
428 fVerbose = flag; 428 fVerbose = flag;
429 } 429 }
430 430
431 G4double G4DNAIRTMoleculeEncounterStepper::Cal 431 G4double G4DNAIRTMoleculeEncounterStepper::CalculateMinTimeStep(G4double currentGlobalTime, G4double definedMinTimeStep){
432 432
433 G4bool start = true; 433 G4bool start = true;
434 G4bool active = false; 434 G4bool active = false;
435 435
436 fUserMinTimeStep = definedMinTimeStep; 436 fUserMinTimeStep = definedMinTimeStep;
437 437
438 if(fReactionSet->Empty()){ 438 if(fReactionSet->Empty()){
439 if(currentGlobalTime == G4Scheduler::I 439 if(currentGlobalTime == G4Scheduler::Instance()->GetStartTime()){
440 440
441 for (auto pTrack : *fpTrackContain 441 for (auto pTrack : *fpTrackContainer->GetMainList())
442 { 442 {
443 if (pTrack == nullptr) 443 if (pTrack == nullptr)
444 { 444 {
445 G4ExceptionDescription exc 445 G4ExceptionDescription exceptionDescription;
446 exceptionDescription << "N 446 exceptionDescription << "No track found.";
447 G4Exception("G4Scheduler:: 447 G4Exception("G4Scheduler::CalculateMinStep", "ITScheduler006",
448 FatalErrorInAr 448 FatalErrorInArgument, exceptionDescription);
449 continue; 449 continue;
450 } 450 }
451 451
452 G4TrackStatus trackStatus = pT 452 G4TrackStatus trackStatus = pTrack->GetTrackStatus();
453 if (trackStatus == fStopAndKil 453 if (trackStatus == fStopAndKill || trackStatus == fStopButAlive)
454 { 454 {
455 start = false; 455 start = false;
456 continue; 456 continue;
457 } 457 }
458 active = true; 458 active = true;
459 } 459 }
460 460
461 if(start){ 461 if(start){
462 return -1; 462 return -1;
463 }if(!active){ 463 }if(!active){
464 G4Scheduler::Instance()->Stop( 464 G4Scheduler::Instance()->Stop();
465 return fSampledMinTimeStep; 465 return fSampledMinTimeStep;
466 } 466 }
467 return fSampledMinTimeStep; 467 return fSampledMinTimeStep;
468 } 468 }
469 for (auto pTrack : *fpTrackContainer-> 469 for (auto pTrack : *fpTrackContainer->GetMainList())
470 { 470 {
471 pTrack->SetGlobalTime(G4Scheduler: 471 pTrack->SetGlobalTime(G4Scheduler::Instance()->GetEndTime());
472 } 472 }
473 return fSampledMinTimeStep; 473 return fSampledMinTimeStep;
474 } 474 }
475 475
476 const auto& fReactionSetInTime = fReaction << 476 auto fReactionSetInTime = fReactionSet->GetReactionsPerTime();
477 fSampledMinTimeStep = fReactionSetInTime.b 477 fSampledMinTimeStep = fReactionSetInTime.begin()->get()->GetTime() - currentGlobalTime;
478 478
479 return fSampledMinTimeStep; 479 return fSampledMinTimeStep;
480 } 480 }
481 481