| 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 // 20/2/2019 25 // 20/2/2019
26 // Author: HoangTRAN 26 // Author: HoangTRAN
27 27
>> 28 #include <memory>
28 #include "G4DNAIndependentReactionTimeStepper. 29 #include "G4DNAIndependentReactionTimeStepper.hh"
29 << 30 #include "G4VDNAReactionModel.hh"
>> 31 #include "G4DNAMolecularReactionTable.hh"
>> 32 #include "G4UnitsTable.hh"
>> 33 #include "G4Molecule.hh"
30 #include "G4ChemicalMoleculeFinder.hh" 34 #include "G4ChemicalMoleculeFinder.hh"
31 #include "G4DNAChemistryManager.hh" 35 #include "G4DNAChemistryManager.hh"
32 #include "G4DNAMakeReaction.hh" 36 #include "G4DNAMakeReaction.hh"
33 #include "G4DNAMolecularReactionTable.hh" <<
34 #include "G4DiffusionControlledReactionModel.h <<
35 #include "G4IRTUtils.hh" <<
36 #include "G4ITReactionChange.hh" 37 #include "G4ITReactionChange.hh"
37 #include "G4Molecule.hh" <<
38 #include "G4Scheduler.hh" 38 #include "G4Scheduler.hh"
39 #include "G4UnitsTable.hh" << 39 #include "G4IRTUtils.hh"
40 #include "G4VDNAReactionModel.hh" <<
41 #include "Randomize.hh" 40 #include "Randomize.hh"
42 << 41 #include "G4DiffusionControlledReactionModel.hh"
43 #include <memory> <<
44 using namespace std; 42 using namespace std;
45 using namespace CLHEP; 43 using namespace CLHEP;
46 44
47 G4DNAIndependentReactionTimeStepper::Utils::Ut << 45 G4DNAIndependentReactionTimeStepper::Utils::Utils(const G4Track& trackA,
48 : fTrackA(trackA), fTrackB(trackB) << 46 const G4Track& trackB)
>> 47 : fTrackA(trackA)
>> 48 , fTrackB(trackB)
49 { 49 {
50 fpMoleculeA = GetMolecule(trackA); 50 fpMoleculeA = GetMolecule(trackA);
51 fpMoleculeB = GetMolecule(trackA); 51 fpMoleculeB = GetMolecule(trackA);
52 } 52 }
53 53
54 G4DNAIndependentReactionTimeStepper::G4DNAInde 54 G4DNAIndependentReactionTimeStepper::G4DNAIndependentReactionTimeStepper()
>> 55 : G4VITTimeStepComputer()
55 { 56 {
56 fReactionSet->SortByTime(); 57 fReactionSet->SortByTime();
57 } 58 }
58 59
59 void G4DNAIndependentReactionTimeStepper::Prep 60 void G4DNAIndependentReactionTimeStepper::Prepare()
60 { 61 {
61 G4VITTimeStepComputer::Prepare(); 62 G4VITTimeStepComputer::Prepare();
62 fSampledPositions.clear(); 63 fSampledPositions.clear();
63 BuildChemicalMoleculeFinder() 64 BuildChemicalMoleculeFinder()
64 } 65 }
65 66
66 void G4DNAIndependentReactionTimeStepper::Init 67 void G4DNAIndependentReactionTimeStepper::InitializeForNewTrack()
67 { 68 {
68 if (fReactants != nullptr) { << 69 if(fReactants != nullptr)
>> 70 {
69 fReactants.reset(); 71 fReactants.reset();
70 } 72 }
71 fSampledMinTimeStep = DBL_MAX; << 73 fSampledMinTimeStep = DBL_MAX;
72 fHasAlreadyReachedNullTime = false; 74 fHasAlreadyReachedNullTime = false;
73 } 75 }
74 76
75 G4double G4DNAIndependentReactionTimeStepper:: << 77 G4double G4DNAIndependentReactionTimeStepper::CalculateStep(
76 << 78 const G4Track& trackA, const G4double& userMinTimeStep)
77 { 79 {
78 auto pMoleculeA = GetMolecule(trackA); 80 auto pMoleculeA = GetMolecule(trackA);
79 InitializeForNewTrack(); 81 InitializeForNewTrack();
80 fUserMinTimeStep = userMinTimeStep; 82 fUserMinTimeStep = userMinTimeStep;
81 fCheckedTracks.insert(trackA.GetTrackID()); 83 fCheckedTracks.insert(trackA.GetTrackID());
82 84
83 #ifdef G4VERBOSE 85 #ifdef G4VERBOSE
84 if (fVerbose != 0) { << 86 if(fVerbose != 0)
>> 87 {
85 G4cout << "_______________________________ 88 G4cout << "________________________________________________________________"
86 "_______" 89 "_______"
87 << G4endl; 90 << G4endl;
88 G4cout << "G4DNAIndependentReactionTimeSte 91 G4cout << "G4DNAIndependentReactionTimeStepper::CalculateStep" << G4endl;
89 G4cout << "Check done for molecule : " << << 92 G4cout << "Check done for molecule : " << pMoleculeA->GetName() << " ("
90 << ") " << G4endl; << 93 << trackA.GetTrackID() << ") " << G4endl;
91 } 94 }
92 #endif 95 #endif
93 96
94 auto pMolConfA = pMoleculeA->GetMolecularCon 97 auto pMolConfA = pMoleculeA->GetMolecularConfiguration();
95 98
96 const auto pReactantList = fMolecularReactio 99 const auto pReactantList = fMolecularReactionTable->CanReactWith(pMolConfA);
97 100
98 if (pReactantList == nullptr) { << 101 if(pReactantList == nullptr)
99 if(fVerbose > 1) { << 102 {
100 G4ExceptionDescription msg; << 103 #ifdef G4VERBOSE
101 msg << "G4DNAIndependentReactionTimeStep << 104 if(fVerbose > 1)
102 "for the reaction because the mol << 105 {
103 << pMoleculeA->GetName() << " does n << 106 G4cout << "!!!!!!!!!!!!!!!!!!!!" << G4endl;
104 << G4endl; << 107 G4cout << "!!! WARNING" << G4endl;
105 G4Exception("G4DNAIndependentReactionTim << 108 G4cout << "G4DNAIndependentReactionTimeStepper::CalculateStep will "
106 "G4DNAIndependentReactionTim << 109 "return infinity "
>> 110 "for the reaction because the molecule "
>> 111 << pMoleculeA->GetName()
>> 112 << " does not have any reactants given in the reaction table."
>> 113 << G4endl;
>> 114 G4cout << "!!!!!!!!!!!!!!!!!!!!" << G4endl;
107 } 115 }
>> 116 #endif
108 return DBL_MAX; 117 return DBL_MAX;
109 } 118 }
110 119
111 auto nbReactives = (G4int)pReactantList->siz << 120 G4int nbReactives = (G4int)pReactantList->size();
112 121
113 if (nbReactives == 0) { << 122 if(nbReactives == 0)
114 if(fVerbose != 0){ << 123 {
115 G4ExceptionDescription msg; << 124 #ifdef G4VERBOSE
116 msg << "G4DNAIndependentReactionTimeStep << 125 // DEBUG
117 "return infinity " << 126 if(fVerbose != 0)
118 "for the reaction because the mol << 127 {
119 << pMoleculeA->GetName() << " does n << 128 G4cout << "!!!!!!!!!!!!!!!!!!!!" << G4endl;
120 << "This message can also result fro << 129 G4cout << "!!! WARNING" << G4endl;
121 "the reaction table." << 130 G4cout << "G4DNAIndependentReactionTimeStepper::CalculateStep will "
122 << G4endl; << 131 "return infinity "
123 G4Exception("G4DNAIndependentReactionTim << 132 "for the reaction because the molecule "
124 "G4DNAIndependentReactionTim << 133 << pMoleculeA->GetName()
>> 134 << " does not have any reactants given in the reaction table."
>> 135 << "This message can also result from a wrong implementation of "
>> 136 "the reaction table."
>> 137 << G4endl;
>> 138 G4cout << "!!!!!!!!!!!!!!!!!!!!" << G4endl;
125 } 139 }
>> 140 #endif
126 return DBL_MAX; 141 return DBL_MAX;
127 } 142 }
128 fReactants = std::make_shared<vector<G4Track 143 fReactants = std::make_shared<vector<G4Track*>>();
129 fReactionModel->Initialise(pMolConfA, trackA 144 fReactionModel->Initialise(pMolConfA, trackA);
130 for (G4int i = 0; i < nbReactives; ++i) { << 145 for(G4int i = 0; i < nbReactives; ++i)
>> 146 {
131 auto pMoleculeB = (*pReactantList)[i]; 147 auto pMoleculeB = (*pReactantList)[i];
132 G4int key = pMoleculeB->GetMoleculeID(); << 148 G4int key = pMoleculeB->GetMoleculeID();
133 149
134 // fRCutOff = G4IRTUtils::GetRCutOff(1 * p 150 // fRCutOff = G4IRTUtils::GetRCutOff(1 * ps);
135 fRCutOff = G4IRTUtils::GetRCutOff(); 151 fRCutOff = G4IRTUtils::GetRCutOff();
136 //________________________________________ 152 //______________________________________________________________
137 // Retrieve reaction range 153 // Retrieve reaction range
138 const G4double Reff = fReactionModel->GetR 154 const G4double Reff = fReactionModel->GetReactionRadius(i);
139 std::vector<std::pair<G4TrackList::iterato 155 std::vector<std::pair<G4TrackList::iterator, G4double>> resultIndices;
140 resultIndices.clear(); 156 resultIndices.clear();
141 G4ChemicalMoleculeFinder::Instance()->Find << 157 G4ChemicalMoleculeFinder::Instance()->FindNearestInRange(
>> 158 trackA, key, fRCutOff, resultIndices);
142 159
143 if (resultIndices.empty()) { << 160 if(resultIndices.empty())
>> 161 {
144 continue; 162 continue;
145 } 163 }
146 for (auto& it : resultIndices) { << 164 for(auto& it : resultIndices)
>> 165 {
147 G4Track* pTrackB = *(std::get<0>(it)); 166 G4Track* pTrackB = *(std::get<0>(it));
148 167
149 if (pTrackB == &trackA) { << 168 if(pTrackB == &trackA)
>> 169 {
150 continue; 170 continue;
151 } 171 }
152 if (pTrackB == nullptr) { << 172 if(pTrackB == nullptr)
>> 173 {
153 G4ExceptionDescription exceptionDescri 174 G4ExceptionDescription exceptionDescription;
154 exceptionDescription << "No trackB no 175 exceptionDescription << "No trackB no valid";
155 G4Exception( << 176 G4Exception("G4DNAIndependentReactionTimeStepper"
156 "G4DNAIndependentReactionTimeStepper << 177 "::CalculateStep()",
157 "::CalculateStep()", << 178 "G4DNAIndependentReactionTimeStepper007", FatalException,
158 "G4DNAIndependentReactionTimeStepper << 179 exceptionDescription);
159 } << 180 }else
160 else { << 181 {
161 if (fCheckedTracks.find(pTrackB->GetTr << 182 if(fCheckedTracks.find(pTrackB->GetTrackID()) != fCheckedTracks.end())
>> 183 {
162 continue; 184 continue;
163 } 185 }
164 186
165 Utils utils(trackA, *pTrackB); 187 Utils utils(trackA, *pTrackB);
166 188
167 auto pMolB = GetMolecule(pTrackB); << 189 auto pMolB = GetMolecule(pTrackB);
168 auto pMolConfB = pMolB->GetMolecularCo << 190 auto pMolConfB = pMolB->GetMolecularConfiguration();
169 G4double distance = (trackA.GetPositio 191 G4double distance = (trackA.GetPosition() - pTrackB->GetPosition()).mag();
170 if (distance * distance < Reff * Reff) << 192 if(distance * distance < Reff * Reff)
171 auto reactionData = fMolecularReacti << 193 {
172 if (G4Scheduler::Instance()->GetGlob << 194 auto reactionData =
173 if (reactionData->GetProbability() << 195 fMolecularReactionTable->GetReactionData(pMolConfA, pMolConfB);
174 if (!fHasAlreadyReachedNullTime) << 196 if(G4Scheduler::Instance()->GetGlobalTime() == G4Scheduler::Instance()->GetStartTime())
>> 197 {
>> 198 if(reactionData->GetProbability() > G4UniformRand())
>> 199 {
>> 200 if(!fHasAlreadyReachedNullTime)
>> 201 {
175 fReactants->clear(); 202 fReactants->clear();
176 fHasAlreadyReachedNullTime = t 203 fHasAlreadyReachedNullTime = true;
177 } 204 }
178 fSampledMinTimeStep = 0.; 205 fSampledMinTimeStep = 0.;
179 CheckAndRecordResults(utils); 206 CheckAndRecordResults(utils);
180 } 207 }
181 } 208 }
182 } 209 }
183 else { << 210 else
>> 211 {
184 G4double tempMinET = GetTimeToEncoun 212 G4double tempMinET = GetTimeToEncounter(trackA, *pTrackB);
185 if (tempMinET < 0 || tempMinET > G4S << 213 if(tempMinET < 0 || tempMinET > G4Scheduler::Instance()->GetEndTime())
>> 214 {
186 continue; 215 continue;
187 } 216 }
188 if (tempMinET >= fSampledMinTimeStep << 217 if(tempMinET >= fSampledMinTimeStep)
>> 218 {
189 continue; 219 continue;
190 } 220 }
191 fSampledMinTimeStep = tempMinET; 221 fSampledMinTimeStep = tempMinET;
192 fReactants->clear(); 222 fReactants->clear();
193 CheckAndRecordResults(utils); 223 CheckAndRecordResults(utils);
194 } 224 }
195 } 225 }
196 } 226 }
197 } 227 }
198 228
199 #ifdef G4VERBOSE 229 #ifdef G4VERBOSE
200 if (fVerbose != 0) { << 230 if(fVerbose != 0)
>> 231 {
201 G4cout << "G4DNAIndependentReactionTimeSte 232 G4cout << "G4DNAIndependentReactionTimeStepper::CalculateStep will finally "
202 "return :" 233 "return :"
203 << G4BestUnit(fSampledMinTimeStep, 234 << G4BestUnit(fSampledMinTimeStep, "Time") << G4endl;
204 235
205 if (fVerbose > 1) { << 236 if(fVerbose > 1)
206 G4cout << "Selected reactants for trackA << 237 {
207 << trackA.GetTrackID() << ") are: << 238 G4cout << "Selected reactants for trackA: " << pMoleculeA->GetName()
>> 239 << " (" << trackA.GetTrackID() << ") are: ";
208 240
209 vector<G4Track*>::iterator it; 241 vector<G4Track*>::iterator it;
210 for (it = fReactants->begin(); it != fRe << 242 for(it = fReactants->begin(); it != fReactants->end(); it++)
>> 243 {
211 G4Track* trackB = *it; 244 G4Track* trackB = *it;
212 G4cout << GetMolecule(trackB)->GetName << 245 G4cout << GetMolecule(trackB)->GetName() << " (" << trackB->GetTrackID()
>> 246 << ") \t ";
213 } 247 }
214 G4cout << G4endl; 248 G4cout << G4endl;
215 } 249 }
216 } 250 }
217 #endif 251 #endif
218 return fSampledMinTimeStep; 252 return fSampledMinTimeStep;
219 } 253 }
220 254
221 void G4DNAIndependentReactionTimeStepper::Chec << 255 void G4DNAIndependentReactionTimeStepper::CheckAndRecordResults(
>> 256 const Utils& utils)
222 { 257 {
223 if (utils.fTrackB.GetTrackStatus() != fAlive << 258 if(utils.fTrackB.GetTrackStatus() != fAlive)
>> 259 {
224 return; 260 return;
225 } 261 }
226 262
227 if (&utils.fTrackB == &utils.fTrackA) { << 263 if(&utils.fTrackB == &utils.fTrackA)
228 G4ExceptionDescription msg; << 264 {
229 msg << "A track is reacting with itself" << 265 G4ExceptionDescription exceptionDescription;
>> 266 exceptionDescription << "A track is reacting with itself"
230 " (which is imposs 267 " (which is impossible) ie fpTrackA == trackB"
231 << G4endl; 268 << G4endl;
232 msg << "Molecule A is of type : " << utils << 269 exceptionDescription << "Molecule A is of type : "
>> 270 << utils.fpMoleculeA->GetName()
233 << " with trackID : " 271 << " with trackID : " << utils.fTrackA.GetTrackID()
234 << " and B : " << uti 272 << " and B : " << utils.fpMoleculeB->GetName()
235 << " with trackID : " << 273 << " with trackID : " << utils.fTrackB.GetTrackID()
>> 274 << G4endl;
236 G4Exception("G4DNAIndependentReactionTimeS 275 G4Exception("G4DNAIndependentReactionTimeStepper::RetrieveResults",
237 "G4DNAIndependentReactionTimeS 276 "G4DNAIndependentReactionTimeStepper003", FatalErrorInArgument,
238 msg); << 277 exceptionDescription);
239 } 278 }
240 279
241 if (fabs(utils.fTrackB.GetGlobalTime() - uti << 280 if(fabs(utils.fTrackB.GetGlobalTime() - utils.fTrackA.GetGlobalTime()) >
242 > utils.fTrackA.GetGlobalTime() * (1. - << 281 utils.fTrackA.GetGlobalTime() * (1. - 1. / 100))
243 { 282 {
244 // DEBUG 283 // DEBUG
245 G4ExceptionDescription msg; << 284 G4ExceptionDescription exceptionDescription;
246 msg << "The interacting tracks are not syn << 285 exceptionDescription
247 msg << "trackB->GetGlobalTime() != fpTrack << 286 << "The interacting tracks are not synchronized in time" << G4endl;
>> 287 exceptionDescription
>> 288 << "trackB->GetGlobalTime() != fpTrackA.GetGlobalTime()" << G4endl;
248 289
249 msg << "fpTrackA : trackID : " << utils.fT << 290 exceptionDescription << "fpTrackA : trackID : "
>> 291 << utils.fTrackA.GetTrackID()
250 << "\t Name :" << uti 292 << "\t Name :" << utils.fpMoleculeA->GetName()
251 << "\t fpTrackA->GetG 293 << "\t fpTrackA->GetGlobalTime() = "
252 << G4BestUnit(utils.f << 294 << G4BestUnit(utils.fTrackA.GetGlobalTime(), "Time")
>> 295 << G4endl;
253 296
254 msg << "trackB : trackID : " << utils.fTra << 297 exceptionDescription << "trackB : trackID : " << utils.fTrackB.GetTrackID()
255 << "\t Name :" << uti 298 << "\t Name :" << utils.fpMoleculeB->GetName()
256 << "\t trackB->GetGlo 299 << "\t trackB->GetGlobalTime() = "
257 << G4BestUnit(utils.f << 300 << G4BestUnit(utils.fTrackB.GetGlobalTime(), "Time")
>> 301 << G4endl;
258 302
259 G4Exception("G4DNAIndependentReactionTimeS 303 G4Exception("G4DNAIndependentReactionTimeStepper::RetrieveResults",
260 "G4DNAIndependentReactionTimeS 304 "G4DNAIndependentReactionTimeStepper004", FatalErrorInArgument,
261 msg); << 305 exceptionDescription);
262 } 306 }
263 fReactants->push_back(const_cast<G4Track*>(& 307 fReactants->push_back(const_cast<G4Track*>(&utils.fTrackB));
264 } 308 }
265 309
266 std::unique_ptr<G4ITReactionChange> G4DNAIndep << 310 std::unique_ptr<G4ITReactionChange>
>> 311 G4DNAIndependentReactionTimeStepper::FindReaction(
267 G4ITReactionSet* pReactionSet, const G4doubl 312 G4ITReactionSet* pReactionSet, const G4double& currentStepTime,
268 const G4double& /*previousStepTime*/, const << 313 const G4double& /*previousStepTime*/,
>> 314 const G4bool& /*reachedUserStepTimeLimit*/)
269 { 315 {
270 if (pReactionSet == nullptr) { << 316 if(pReactionSet == nullptr)
>> 317 {
271 return nullptr; 318 return nullptr;
272 } 319 }
273 320
274 G4ITReactionPerTime& reactionPerTime = pReac 321 G4ITReactionPerTime& reactionPerTime = pReactionSet->GetReactionsPerTime();
275 if (reactionPerTime.empty()) { << 322 if(reactionPerTime.empty())
>> 323 {
276 return nullptr; 324 return nullptr;
277 } 325 }
278 for (auto reaction_i = reactionPerTime.begin <<
279 reaction_i = reactionPerTime.begin()) <<
280 { <<
281 if ((*reaction_i)->GetTime() > currentStep <<
282 fReactionSet->CleanAllReaction(); <<
283 return nullptr; <<
284 } <<
285 326
>> 327 for(auto reaction_i = reactionPerTime.begin();
>> 328 reaction_i != reactionPerTime.end(); reaction_i = reactionPerTime.begin())
>> 329 {
286 G4Track* pTrackA = (*reaction_i)->GetReact 330 G4Track* pTrackA = (*reaction_i)->GetReactants().first;
287 if (pTrackA->GetTrackStatus() == fStopAndK << 331 if(pTrackA->GetTrackStatus() == fStopAndKill)
>> 332 {
288 continue; 333 continue;
289 } 334 }
290 G4Track* pTrackB = (*reaction_i)->GetReact 335 G4Track* pTrackB = (*reaction_i)->GetReactant(pTrackA);
291 if (pTrackB->GetTrackStatus() == fStopAndK << 336 if(pTrackB->GetTrackStatus() == fStopAndKill)
>> 337 {
292 continue; 338 continue;
293 } 339 }
294 340
295 if (pTrackB == pTrackA) { << 341 if(pTrackB == pTrackA)
296 G4ExceptionDescription msg; << 342 {
297 msg << "The IT reaction process sent bac << 343 G4ExceptionDescription exceptionDescription;
>> 344 exceptionDescription << "The IT reaction process sent back a reaction "
298 "between trackA 345 "between trackA and trackB. ";
299 msg << "The problem is trackA == trackB" << 346 exceptionDescription << "The problem is trackA == trackB";
300 G4Exception("G4DNAIndependentReactionTim 347 G4Exception("G4DNAIndependentReactionTimeStepper::FindReaction",
301 "G4DNAIndependentReactionTim 348 "G4DNAIndependentReactionTimeStepper02", FatalErrorInArgument,
302 msg); << 349 exceptionDescription);
303 } 350 }
304 pReactionSet->SelectThisReaction(*reaction 351 pReactionSet->SelectThisReaction(*reaction_i);
305 if (fpReactionProcess != nullptr << 352 if(fpReactionProcess != nullptr &&
306 && fpReactionProcess->TestReactibility << 353 fpReactionProcess->TestReactibility(*pTrackA, *pTrackB, currentStepTime,
>> 354 false))
307 { 355 {
308 if ((fSampledPositions.find(pTrackA->Get << 356 if((fSampledPositions.find(pTrackA->GetTrackID()) ==
309 && (fSampledPositions.find(pTrackB- << 357 fSampledPositions.end() &&
>> 358 (fSampledPositions.find(pTrackB->GetTrackID()) ==
>> 359 fSampledPositions.end())))
310 { 360 {
311 G4ExceptionDescription msg; << 361 G4ExceptionDescription exceptionDescription;
312 msg << "The positions of trackA and tr << 362 exceptionDescription
>> 363 << "The positions of trackA and trackB have no counted ";
313 G4Exception("G4DNAIndependentReactionT 364 G4Exception("G4DNAIndependentReactionTimeStepper::FindReaction",
314 "G4DNAIndependentReactionT << 365 "G4DNAIndependentReactionTimeStepper0001",
315 msg); << 366 FatalErrorInArgument, exceptionDescription);
316 } 367 }
317 368
318 pTrackA->SetPosition(fSampledPositions[p 369 pTrackA->SetPosition(fSampledPositions[pTrackA->GetTrackID()]);
319 pTrackB->SetPosition(fSampledPositions[p 370 pTrackB->SetPosition(fSampledPositions[pTrackB->GetTrackID()]);
320 auto pReactionChange = fpReactionProcess << 371 auto pReactionChange =
321 if (pReactionChange == nullptr) { << 372 fpReactionProcess->MakeReaction(*pTrackA, *pTrackB);
>> 373 if(pReactionChange == nullptr)
>> 374 {
322 return nullptr; 375 return nullptr;
323 } 376 }
324 return pReactionChange; 377 return pReactionChange;
325 } 378 }
326 } 379 }
327 return nullptr; 380 return nullptr;
328 } 381 }
329 382
330 void G4DNAIndependentReactionTimeStepper::SetR << 383 void G4DNAIndependentReactionTimeStepper::SetReactionModel(
>> 384 G4VDNAReactionModel* pReactionModel)
331 { 385 {
332 fReactionModel = pReactionModel; 386 fReactionModel = pReactionModel;
333 } 387 }
334 388
335 G4VDNAReactionModel* G4DNAIndependentReactionT 389 G4VDNAReactionModel* G4DNAIndependentReactionTimeStepper::GetReactionModel()
336 { 390 {
337 return fReactionModel; 391 return fReactionModel;
338 } 392 }
339 393
340 void G4DNAIndependentReactionTimeStepper::SetV 394 void G4DNAIndependentReactionTimeStepper::SetVerbose(G4int flag)
341 { 395 {
342 fVerbose = flag; 396 fVerbose = flag;
343 } 397 }
344 398
345 G4double G4DNAIndependentReactionTimeStepper:: << 399 G4double G4DNAIndependentReactionTimeStepper::GetTimeToEncounter(
346 << 400 const G4Track& trackA, const G4Track& trackB)
347 { 401 {
348 G4double timeToReaction = dynamic_cast<G4Dif << 402 G4double timeToReaction =
349 ->GetTimeToEncou << 403 dynamic_cast<G4DiffusionControlledReactionModel*>(fReactionModel)
>> 404 ->GetTimeToEncounter(trackA, trackB);
350 return timeToReaction; 405 return timeToReaction;
351 } 406 }
352 407
353 void G4DNAIndependentReactionTimeStepper::SetR << 408 void G4DNAIndependentReactionTimeStepper::SetReactionProcess(
>> 409 G4VITReactionProcess* pReactionProcess)
354 { 410 {
355 fpReactionProcess = pReactionProcess; 411 fpReactionProcess = pReactionProcess;
356 } 412 }
357 G4double G4DNAIndependentReactionTimeStepper:: << 413 G4double G4DNAIndependentReactionTimeStepper::CalculateMinTimeStep(
358 << 414 G4double /*currentGlobalTime*/, G4double definedMinTimeStep)
359 { 415 {
360 G4double fTSTimeStep = DBL_MAX; 416 G4double fTSTimeStep = DBL_MAX;
361 fCheckedTracks.clear(); 417 fCheckedTracks.clear();
362 418
363 for (auto pTrack : *fpTrackContainer->GetMai << 419 for(auto pTrack : *fpTrackContainer->GetMainList())
364 if (pTrack == nullptr) { << 420 {
365 G4ExceptionDescription msg; << 421 if(pTrack == nullptr)
366 msg << "No track found."; << 422 {
>> 423 G4ExceptionDescription exceptionDescription;
>> 424 exceptionDescription << "No track found.";
367 G4Exception("G4DNAIndependentReactionTim 425 G4Exception("G4DNAIndependentReactionTimeStepper::CalculateMinTimeStep",
368 "G4DNAIndependentReactionTim << 426 "G4DNAIndependentReactionTimeStepper006",
369 msg); << 427 FatalErrorInArgument, exceptionDescription);
370 continue; 428 continue;
371 } 429 }
372 430
373 G4TrackStatus trackStatus = pTrack->GetTra 431 G4TrackStatus trackStatus = pTrack->GetTrackStatus();
374 if (trackStatus == fStopAndKill || trackSt << 432 if(trackStatus == fStopAndKill || trackStatus == fStopButAlive)
>> 433 {
375 continue; 434 continue;
376 } 435 }
377 436
378 G4double sampledMinTimeStep = CalculateSte << 437 G4double sampledMinTimeStep = CalculateStep(*pTrack, definedMinTimeStep);
379 G4TrackVectorHandle reactants = GetReactan 438 G4TrackVectorHandle reactants = GetReactants();
380 439
381 if (sampledMinTimeStep < fTSTimeStep) { << 440 if(sampledMinTimeStep < fTSTimeStep)
>> 441 {
382 fTSTimeStep = sampledMinTimeStep; 442 fTSTimeStep = sampledMinTimeStep;
383 if (reactants) { << 443 fReactionSet->CleanAllReaction();
384 fReactionSet->AddReactions(fTSTimeStep << 444 if(reactants)
>> 445 {
>> 446 fReactionSet->AddReactions(fTSTimeStep, const_cast<G4Track*>(pTrack),
>> 447 reactants);
385 448
386 fSampledPositions[pTrack->GetTrackID() 449 fSampledPositions[pTrack->GetTrackID()] = pTrack->GetPosition();
387 for (const auto& it : *fReactants) { << 450 for(const auto& it : *fReactants)
>> 451 {
388 auto pTrackB = it; 452 auto pTrackB = it;
>> 453 // G4cout<<"position : "<<pTrackB->GetTrackID()<<G4endl;
389 fSampledPositions[pTrackB->GetTrackI 454 fSampledPositions[pTrackB->GetTrackID()] = pTrackB->GetPosition();
390 } 455 }
391 ResetReactants(); 456 ResetReactants();
392 } 457 }
393 } 458 }
394 else if (fTSTimeStep == sampledMinTimeStep << 459 else if(fTSTimeStep == sampledMinTimeStep && G4bool(reactants))
395 fReactionSet->AddReactions(fTSTimeStep, << 460 {
>> 461 fReactionSet->AddReactions(fTSTimeStep, const_cast<G4Track*>(pTrack),
>> 462 reactants);
396 463
397 fSampledPositions[pTrack->GetTrackID()] 464 fSampledPositions[pTrack->GetTrackID()] = pTrack->GetPosition();
398 for (const auto& it : *fReactants) { << 465 for(const auto& it : *fReactants)
>> 466 {
399 auto pTrackB = it; 467 auto pTrackB = it;
>> 468 // G4cout<<"position : "<<pTrackB->GetTrackID()<<G4endl;
400 fSampledPositions[pTrackB->GetTrackID( 469 fSampledPositions[pTrackB->GetTrackID()] = pTrackB->GetPosition();
401 } 470 }
402 ResetReactants(); 471 ResetReactants();
403 } 472 }
404 else if (reactants) { << 473 else if(reactants)
>> 474 {
405 ResetReactants(); 475 ResetReactants();
406 } 476 }
407 } 477 }
>> 478
408 return fTSTimeStep; 479 return fTSTimeStep;
409 } 480 }
410 481