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