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