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