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