| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 //
27 // Author: Mathieu Karamitros (kara@cenbg.in2p
28 //
29 // WARNING : This class is released as a proto
30 // It might strongly evolve or even disapear i
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(co
61 co
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
69 }
70
71 G4DNAMoleculeEncounterStepper::G4DNAMoleculeEn
72 :
73 fMolecularReactionTable(reference_cast<co
74 {
75 fpTrackContainer = G4ITTrackHolder::Instan
76 fReactionSet = G4ITReactionSet::Instance()
77 }
78
79 G4DNAMoleculeEncounterStepper::~G4DNAMoleculeE
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()->UpdatePositi
93
94 #if defined (DEBUG_MEM)
95 mem_second = MemoryUsage();
96 mem_diff = mem_second - mem_first;
97 G4cout << "\t || MEM || G4DNAMoleculeEncou
98 "After computing G4ITManager<G4Molecul
99 "UpdatePositionMap, diff is : " << mem
100 #endif
101 }
102
103 void G4DNAMoleculeEncounterStepper::Initialize
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_infinit
117 && value == std::numeric_limits<T>::in
118 }
119
120 G4double
121 G4DNAMoleculeEncounterStepper::CalculateStep(c
122 c
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 << "G4DNAMoleculeEncounterStepp
135 G4cout << "Check done for molecule : "
136 << " (" << trackA.GetTrackID() <<
137 << G4endl;
138 }
139 #endif
140
141 //________________________________________
142 // Retrieve general informations for makin
143 auto pMolConfA = pMoleculeA->GetMolecularC
144
145 const auto pReactantList = fMolecularReact
146
147 if (pReactantList == nullptr)
148 {
149 #ifdef G4VERBOSE
150 // DEBUG
151 if (fVerbose > 1)
152 {
153 G4cout << "!!!!!!!!!!!!!!!!!!!!" <
154 G4cout << "!!! WARNING" << G4endl;
155 G4cout << "G4MoleculeEncounterStep
156 "for the reaction because the
157 << pMoleculeA->GetName()
158 << " does not have any reactan
159 << G4endl;
160 G4cout << "!!!!!!!!!!!!!!!!!!!!" <
161 }
162 #endif
163 return DBL_MAX;
164 }
165
166 auto nbReactives = (G4int)pReactantList->
167
168 if (nbReactives == 0)
169 {
170 #ifdef G4VERBOSE
171 // DEBUG
172 if (fVerbose != 0)
173 {
174 // TODO replace with the warning m
175 G4cout << "!!!!!!!!!!!!!!!!!!!!" <
176 G4cout << "!!! WARNING" << G4endl;
177 G4cout << "G4MoleculeEncounterStep
178 "for the reaction because the
179 << pMoleculeA->GetName()
180 << " does not have any reactan
181 << "This message can also resu
182 << G4endl;
183 G4cout << "!!!!!!!!!!!!!!!!!!!!" <
184 }
185 #endif
186 return DBL_MAX;
187 }
188
189 fReactants = std::make_shared<vector<G4Tra
190 fReactionModel->Initialise(pMolConfA, trac
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->Get
201
202 //____________________________________
203 // Use KdTree algorithm to find closes
204 G4KDTreeResultHandle resultsNearest(
205 G4MoleculeFinder::Instance()->Find
206
207
208 if (static_cast<int>(resultsNearest) =
209
210 G4double r2 = resultsNearest->GetDista
211 Utils utils(trackA, pMoleculeB);
212
213 if (r2 <= R * R) // ==> Record in rang
214 {
215 // Entering in this condition may
216 // to each other
217 // Therefore, if we only take the
218 // reacted. Instead, we will take
219
220 if (!fHasAlreadyReachedNullTime)
221 {
222 fReactants->clear();
223 fHasAlreadyReachedNullTime = t
224 }
225
226 fSampledMinTimeStep = 0.;
227 G4KDTreeResultHandle resultsInRang
228 G4MoleculeFinder::Instance()->
229
230
231 CheckAndRecordResults(utils,
232 #ifdef G4VERBOSE
233 R,
234 #endif
235 resultsInRan
236 }
237 else
238 {
239 G4double r = sqrt(r2);
240 G4double tempMinET = pow(r - R, 2)
241 // constant = 16 * (fDA + fDB + 2*
242
243 if (tempMinET <= fSampledMinTimeSt
244 {
245 if (fUserMinTimeStep < DBL_MAX
246 && tempMinET <= fUserMinTi
247 {
248 if (fSampledMinTimeStep >
249 {
250 fReactants->clear();
251 }
252
253 fSampledMinTimeStep = fUse
254
255 G4double range = R + sqrt(
256
257 G4KDTreeResultHandle resul
258 G4MoleculeFinder::Inst
259 FindNearestInRange(pMo
260 pMo
261 ran
262
263 CheckAndRecordResults(util
264 #ifdef G4VERBOSE
265 rang
266 #endif
267 resu
268 }
269 else // ==> Record nearest
270 {
271 if (tempMinET < fSampledMi
272 // to avoid cases wher
273 {
274 fSampledMinTimeStep =
275 fReactants->clear();
276 }
277
278 CheckAndRecordResults(util
279 #ifdef G4VERBOSE
280 R,
281 #endif
282 resu
283 }
284 }
285 }
286 }
287
288 #ifdef G4VERBOSE
289 if (fVerbose != 0)
290 {
291 G4cout << "G4MoleculeEncounterStepper:
292 << G4BestUnit(fSampledMinTimeStep,
293
294 if (fVerbose > 1)
295 {
296 G4cout << "Selected reactants for
297 << " (" << trackA.GetTrackID()
298
299 vector<G4Track*>::iterator it;
300 for (it = fReactants->begin(); it
301 {
302 G4Track* trackB = *it;
303 G4cout << GetMolecule(trackB)-
304 << trackB->GetTrackID() <<
305 }
306 G4cout << G4endl;
307 }
308 }
309 #endif
310 return fSampledMinTimeStep;
311 }
312
313 void G4DNAMoleculeEncounterStepper::CheckAndRe
314 #ifdef G4VERBOSE
315
316 #endif
317
318 {
319 if (static_cast<int>(results) == 0)
320 {
321 #ifdef G4VERBOSE
322 if (fVerbose > 1)
323 {
324 G4cout << "No molecule " << utils.
325 << " found to react with " <<
326 << G4endl;
327 }
328 #endif
329 return;
330 }
331
332 for (results->Rewind(); !results->End(); r
333 {
334 auto reactiveB = results->GetItem<G4I
335
336 if (reactiveB == nullptr)
337 {
338 continue;
339 }
340
341 G4Track *trackB = reactiveB->GetTrack(
342
343 if (trackB == nullptr)
344 {
345 G4ExceptionDescription exceptionDe
346 exceptionDescription
347 << "The reactant B found using
348 "track attached to it. If this
349 "not record this molecule in t
350 << G4endl;
351 G4Exception("G4DNAMoleculeEncounte
352 "MoleculeEncounterStep
353 exceptionDescription);
354 continue;
355 }
356
357 if (trackB->GetTrackStatus() != fAlive
358 {
359 continue;
360 }
361
362 if (trackB == &utils.fpTrackA)
363 {
364 G4ExceptionDescription exceptionDe
365 exceptionDescription
366 << "A track is reacting with i
367 << G4endl;
368 exceptionDescription << "Molecule
369 << utils.fpMoleculeA->GetName(
370 << utils.fpTrackA.GetTrackID()
371
372 G4Exception("G4DNAMoleculeEncounte
373 "MoleculeEncounterStep
374 exceptionDescription);
375
376 }
377
378 if (fabs(trackB->GetGlobalTime() - uti
379 > utils.fpTrackA.GetGlobalTime() * (1
380 {
381 // DEBUG
382 G4ExceptionDescription exceptionDe
383 exceptionDescription
384 << "The interacting tracks are
385 exceptionDescription
386 << "trackB->GetGlobalTime() !=
387
388 exceptionDescription << "fpTrackA
389 << "\t Name :" << utils.fpMole
390 << "\t fpTrackA->GetGlobalTime
391 << G4BestUnit(utils.fpTrackA.G
392
393 exceptionDescription << "trackB :
394 << "\t Name :" << utils.fpMole
395 << "\t trackB->GetGlobalTime()
396 << G4BestUnit(trackB->GetGloba
397
398 G4Exception("G4DNAMoleculeEncounte
399 "MoleculeEncounterStep
400 exceptionDescription);
401 }
402
403 #ifdef G4VERBOSE
404 if (fVerbose > 1)
405 {
406
407 G4double r2 = results->GetDistance
408 G4cout << "\t ********************
409 G4cout << "\t Reaction between "
410 << utils.fpMoleculeA->GetName(
411 << " & " << utils.fpMoleculeB-
412 << "Interaction Range = "
413 << G4BestUnit(R, "Length") <<
414 G4cout << "\t Real distance betwee
415 << G4BestUnit((utils.fpTrackA.
416 G4cout << "\t Distance between rea
417 << G4BestUnit(sqrt(r2), "Lengt
418
419 }
420 #endif
421
422 fReactants->push_back(trackB);
423 }
424 }
425
426 void G4DNAMoleculeEncounterStepper::SetReactio
427 {
428 fReactionModel = pReactionModel;
429 }
430
431 G4VDNAReactionModel* G4DNAMoleculeEncounterSte
432 {
433 return fReactionModel;
434 }
435
436 void G4DNAMoleculeEncounterStepper::SetVerbose
437 {
438 fVerbose = flag;
439 }
440
441 G4double G4DNAMoleculeEncounterStepper::Calcul
442
443 G4double fTSTimeStep = DBL_MAX;
444
445 for (auto pTrack : *fpTrackContainer->GetM
446 {
447 if (pTrack == nullptr)
448 {
449 G4ExceptionDescription exceptionDe
450 exceptionDescription << "No track
451 G4Exception("G4Scheduler::Calculat
452 FatalErrorInArgument,
453 continue;
454 }
455
456 G4TrackStatus trackStatus = pTrack->Ge
457 if (trackStatus == fStopAndKill || tra
458 {
459 continue;
460 }
461
462 G4double sampledMinTimeStep = Calculat
463 G4TrackVectorHandle reactants = GetRea
464
465 if (sampledMinTimeStep < fTSTimeStep)
466 {
467 fTSTimeStep = sampledMinTimeStep;
468 fReactionSet->CleanAllReaction();
469 if (reactants)
470 {
471 fReactionSet->AddReactions(fTS
472 con
473 std
474 ResetReactants();
475 }
476 }
477 else if (fTSTimeStep == sampledMinTim
478 {
479 fReactionSet->AddReactions(fTSTim
480 const_
481 std::m
482 ResetReactants();
483 }
484 else if (reactants)
485 {
486 ResetReactants();
487 }
488 }
489
490 return fTSTimeStep;
491 }
492