| 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 (AT) cenbg 27 // Author: Mathieu Karamitros (kara (AT) cenbg . in2p3 . fr)
28 // 28 //
29 // History: 29 // History:
30 // ----------- 30 // -----------
31 // 10 Oct 2011 M.Karamitros created 31 // 10 Oct 2011 M.Karamitros created
32 // 32 //
33 // ------------------------------------------- 33 // -------------------------------------------------------------------
34 34
35 #include "G4ITModelProcessor.hh" 35 #include "G4ITModelProcessor.hh"
36 #include "G4VITTimeStepComputer.hh" 36 #include "G4VITTimeStepComputer.hh"
37 #include "G4VITReactionProcess.hh" 37 #include "G4VITReactionProcess.hh"
38 #include "G4ITReaction.hh" 38 #include "G4ITReaction.hh"
39 #include "G4ITTrackHolder.hh" 39 #include "G4ITTrackHolder.hh"
40 #include "G4ITTrackingManager.hh" 40 #include "G4ITTrackingManager.hh"
41 #include "G4VITStepModel.hh" 41 #include "G4VITStepModel.hh"
42 #include "G4UserTimeStepAction.hh" 42 #include "G4UserTimeStepAction.hh"
43 #include "G4UnitsTable.hh" 43 #include "G4UnitsTable.hh"
44 #include "G4Scheduler.hh" <<
45 #include "G4SystemOfUnits.hh" <<
46 #include <vector> 44 #include <vector>
47 45
48 //#define DEBUG_MEM 46 //#define DEBUG_MEM
49 47
50 #ifdef DEBUG_MEM 48 #ifdef DEBUG_MEM
51 #include "G4MemStat.hh" 49 #include "G4MemStat.hh"
52 using namespace G4MemStat; 50 using namespace G4MemStat;
53 #endif 51 #endif
54 52
55 G4ITModelProcessor::G4ITModelProcessor() 53 G4ITModelProcessor::G4ITModelProcessor()
56 { 54 {
57 fpTrack = nullptr; 55 fpTrack = nullptr;
58 fInitialized = false; 56 fInitialized = false;
59 fUserMinTimeStep = -1.; 57 fUserMinTimeStep = -1.;
60 fTSTimeStep = DBL_MAX; 58 fTSTimeStep = DBL_MAX;
61 fpTrackingManager = nullptr; 59 fpTrackingManager = nullptr;
62 fReactionSet = nullptr; 60 fReactionSet = nullptr;
63 fpTrackContainer = nullptr; 61 fpTrackContainer = nullptr;
64 fpModelHandler = nullptr; 62 fpModelHandler = nullptr;
65 fpActiveModelWithMinTimeStep = nullptr; 63 fpActiveModelWithMinTimeStep = nullptr;
66 fComputeTimeStep = false; 64 fComputeTimeStep = false;
67 fComputeReaction = false; 65 fComputeReaction = false;
68 } 66 }
69 67
70 G4ITModelProcessor::~G4ITModelProcessor() = de 68 G4ITModelProcessor::~G4ITModelProcessor() = default;
71 69
72 void G4ITModelProcessor::RegisterModel(double 70 void G4ITModelProcessor::RegisterModel(double time, G4VITStepModel* model)
73 { 71 {
74 fpModelHandler->RegisterModel(model, time) 72 fpModelHandler->RegisterModel(model, time);
75 } 73 }
76 74
77 void G4ITModelProcessor::Initialize() 75 void G4ITModelProcessor::Initialize()
78 { 76 {
79 fpModelHandler->Initialize(); 77 fpModelHandler->Initialize();
80 fReactionSet = G4ITReactionSet::Instance() 78 fReactionSet = G4ITReactionSet::Instance();
81 fpTrackContainer = G4ITTrackHolder::Instan 79 fpTrackContainer = G4ITTrackHolder::Instance();
82 fInitialized = true; 80 fInitialized = true;
83 fComputeTimeStep = false; 81 fComputeTimeStep = false;
84 fComputeReaction = false; 82 fComputeReaction = false;
85 if (fpModelHandler->GetTimeStepComputerFla 83 if (fpModelHandler->GetTimeStepComputerFlag())
86 { 84 {
87 fComputeTimeStep = true; 85 fComputeTimeStep = true;
88 } 86 }
89 if (fpModelHandler->GetReactionProcessFlag 87 if (fpModelHandler->GetReactionProcessFlag())
90 { 88 {
91 fComputeReaction = true; 89 fComputeReaction = true;
92 } 90 }
93 } 91 }
94 92
95 G4double G4ITModelProcessor::CalculateMinTimeS 93 G4double G4ITModelProcessor::CalculateMinTimeStep(G4double currentGlobalTime,
96 94 G4double definedMinTimeStep)
97 { 95 {
98 96
99 #if defined (DEBUG_MEM) && defined (DEBUG_MEM_ 97 #if defined (DEBUG_MEM) && defined (DEBUG_MEM_DETAILED_STEPPING)
100 MemStat mem_first, mem_second, mem_diff; 98 MemStat mem_first, mem_second, mem_diff;
101 mem_first = MemoryUsage(); 99 mem_first = MemoryUsage();
102 #endif 100 #endif
103 101
104 fpActiveModelWithMinTimeStep = nullptr; 102 fpActiveModelWithMinTimeStep = nullptr;
105 fTSTimeStep = DBL_MAX; 103 fTSTimeStep = DBL_MAX;
106 104
107 InitializeStepper(currentGlobalTime, defin 105 InitializeStepper(currentGlobalTime, definedMinTimeStep);
108 106
109 #if defined (DEBUG_MEM) && defined (DEBUG_MEM_ 107 #if defined (DEBUG_MEM) && defined (DEBUG_MEM_DETAILED_STEPPING)
110 mem_second = MemoryUsage(); 108 mem_second = MemoryUsage();
111 mem_diff = mem_second-mem_first; 109 mem_diff = mem_second-mem_first;
112 G4cout << "\t || MEM || G4Scheduler::Calcu 110 G4cout << "\t || MEM || G4Scheduler::CalculateMinTimeStep || After "
113 "computing fpModelProcessor -> InitializeS 111 "computing fpModelProcessor -> InitializeStepper, diff is : "
114 << mem_diff 112 << mem_diff
115 << G4endl; 113 << G4endl;
116 #endif 114 #endif
117 115
118 for (auto& pStepModel : fActiveModels) << 116 for (auto pTrack : *fpTrackContainer->GetMainList())
119 { 117 {
120 fTSTimeStep = << 118 if (pTrack == nullptr)
121 pStepModel->GetTimeStepper()->Calc << 119 {
122 currentGlobalTime, << 120 G4ExceptionDescription exceptionDescription;
123 definedMinTimeStep); << 121 exceptionDescription << "No track found.";
124 << 122 G4Exception("G4Scheduler::CalculateMinStep", "ITScheduler006",
125 fpActiveModelWithMinTimeStep = pStepMo << 123 FatalErrorInArgument, exceptionDescription);
126 << 124 continue;
127 if(fTSTimeStep == -1){ <<
128 fpActiveModelWithMinTimeStep->GetR <<
129 if(fReactionSet->Empty()) return D <<
130 const auto& fReactionSetInTime = f <<
131 fTSTimeStep = fReactionSetInTime.b <<
132 } 125 }
>> 126
>> 127 #ifdef DEBUG
>> 128 G4cout << "*_* " << GetIT(track)->GetName()
>> 129 << " ID: " << track->GetTrackID()
>> 130 << " at time : " << track->GetGlobalTime()
>> 131 << G4endl;
>> 132 #endif
>> 133
>> 134 G4TrackStatus trackStatus = pTrack->GetTrackStatus();
>> 135 if (trackStatus == fStopAndKill || trackStatus == fStopButAlive)
>> 136 {
>> 137 continue;
>> 138 }
>> 139
>> 140 CalculateTimeStep(pTrack, definedMinTimeStep);
>> 141 // if MT mode at track level, this command should be displaced
>> 142 ExtractTimeStepperData();
133 } 143 }
134 144
135 #if defined (DEBUG_MEM) && defined (DEBUG_MEM_ 145 #if defined (DEBUG_MEM) && defined (DEBUG_MEM_DETAILED_STEPPING)
136 mem_second = MemoryUsage(); 146 mem_second = MemoryUsage();
137 mem_diff = mem_second-mem_first; 147 mem_diff = mem_second-mem_first;
138 G4cout << "\t || MEM || G4Scheduler::Calcu 148 G4cout << "\t || MEM || G4Scheduler::CalculateMinTimeStep || "
139 "After looping on tracks, diff is : " << m 149 "After looping on tracks, diff is : " << mem_diff << G4endl;
140 #endif 150 #endif
141 return fTSTimeStep; 151 return fTSTimeStep;
142 } 152 }
143 153
>> 154 //_________________________________________________________________________
>> 155
>> 156 void G4ITModelProcessor::ExtractTimeStepperData()
>> 157 {
>> 158 if (fpTrack == nullptr)
>> 159 {
>> 160 CleanProcessor();
>> 161 return;
>> 162 }
>> 163
>> 164 for (auto pStepModel : fActiveModels)
>> 165 {
>> 166 if (pStepModel == nullptr)
>> 167 {
>> 168 continue;
>> 169 }
>> 170
>> 171 auto pTimeStepper = pStepModel->GetTimeStepper();
>> 172 G4double sampledMinTimeStep = pTimeStepper->GetSampledMinTimeStep();
>> 173 G4TrackVectorHandle reactants = pTimeStepper->GetReactants();
>> 174
>> 175 if (sampledMinTimeStep < fTSTimeStep)
>> 176 {
>> 177 fpActiveModelWithMinTimeStep = pStepModel;
>> 178 fTSTimeStep = sampledMinTimeStep;
>> 179 //fReactingTracks.clear();
>> 180
>> 181 fReactionSet->CleanAllReaction();
>> 182 if (reactants)
>> 183 {
>> 184 // fReactingTracks.insert(make_pair(track, reactants));
>> 185 fReactionSet->AddReactions(fTSTimeStep,
>> 186 const_cast<G4Track*>(fpTrack),
>> 187 reactants);
>> 188 pTimeStepper->ResetReactants();
>> 189 }
>> 190 }
>> 191 else if (fTSTimeStep == sampledMinTimeStep && bool(reactants))
>> 192 {
>> 193 // fReactingTracks.insert(make_pair(track, reactants));
>> 194 fReactionSet->AddReactions(fTSTimeStep,
>> 195 const_cast<G4Track*>(fpTrack),
>> 196 reactants);
>> 197 pTimeStepper->ResetReactants();
>> 198 }
>> 199 else if (reactants)
>> 200 {
>> 201 pTimeStepper->ResetReactants();
>> 202 }
>> 203 }
>> 204
>> 205 CleanProcessor();
>> 206 }
>> 207
144 //____________________________________________ 208 //______________________________________________________________________________
145 209
146 void G4ITModelProcessor::InitializeStepper(G4d 210 void G4ITModelProcessor::InitializeStepper(G4double currentGlobalTime,
147 G4d 211 G4double userMinTime)
148 { 212 {
149 G4VITTimeStepComputer::SetTimes(currentGlo 213 G4VITTimeStepComputer::SetTimes(currentGlobalTime, userMinTime);
150 214
151 #if defined (DEBUG_MEM) 215 #if defined (DEBUG_MEM)
152 MemStat mem_first, mem_second, mem_diff; 216 MemStat mem_first, mem_second, mem_diff;
153 mem_first = MemoryUsage(); 217 mem_first = MemoryUsage();
154 #endif 218 #endif
155 219
156 fActiveModels = fpModelHandler->GetActiveM 220 fActiveModels = fpModelHandler->GetActiveModels(currentGlobalTime);
157 221
158 for (auto& pModel : fActiveModels) 222 for (auto& pModel : fActiveModels)
159 { 223 {
160 pModel->PrepareNewTimeStep(); 224 pModel->PrepareNewTimeStep();
161 } 225 }
162 226
163 #if defined (DEBUG_MEM) 227 #if defined (DEBUG_MEM)
164 mem_second = MemoryUsage(); 228 mem_second = MemoryUsage();
165 mem_diff = mem_second-mem_first; 229 mem_diff = mem_second-mem_first;
166 G4cout << "\t || MEM || G4ITModelP 230 G4cout << "\t || MEM || G4ITModelProcessor::InitializeStepper || After computing stepper -> Prepare(), diff is : " << mem_diff << G4endl;
167 #endif 231 #endif
168 232
169 } 233 }
170 234
>> 235 //______________________________________________________________________________
>> 236 void G4ITModelProcessor::CalculateTimeStep(const G4Track* pTrack,
>> 237 const G4double userMinTimeStep)
>> 238 {
>> 239 CleanProcessor();
>> 240 if (pTrack == nullptr)
>> 241 {
>> 242 G4ExceptionDescription exceptionDescription;
>> 243 exceptionDescription << "No track was passed to the method.";
>> 244 G4Exception("G4ITModelProcessor::CalculateStep",
>> 245 "ITModelProcessor004",
>> 246 FatalErrorInArgument,
>> 247 exceptionDescription);
>> 248 }
>> 249 SetTrack(pTrack);
>> 250 fUserMinTimeStep = userMinTimeStep;
>> 251
>> 252 DoCalculateStep();
>> 253 }
>> 254
>> 255 //______________________________________________________________________________
>> 256
>> 257 void G4ITModelProcessor::DoCalculateStep()
>> 258 {
>> 259 for (auto& pStepModel : fActiveModels)
>> 260 {
>> 261 pStepModel->GetTimeStepper()->CalculateStep(*fpTrack, fUserMinTimeStep);
>> 262 }
>> 263 }
>> 264
171 //____________________________________________ 265 //_________________________________________________________________________
172 266
173 void G4ITModelProcessor::ComputeTrackReaction( 267 void G4ITModelProcessor::ComputeTrackReaction(G4ITStepStatus fITStepStatus,
174 268 G4double fGlobalTime,
175 269 G4double currentTimeStep,
176 << 270 G4double previousTimeStep,
177 271 G4bool reachedUserTimeLimit,
178 272 G4double fTimeTolerance,
179 273 G4UserTimeStepAction* fpUserTimeStepAction,
180 274 G4int
181 #ifdef G4VERBOSE 275 #ifdef G4VERBOSE
182 fVerbose 276 fVerbose
183 #endif 277 #endif
184 ) 278 )
185 { 279 {
>> 280 // if (fReactingTracks.empty())
186 if (fReactionSet->Empty()) 281 if (fReactionSet->Empty())
187 { 282 {
188 return; 283 return;
189 } 284 }
190 285
191 if (fITStepStatus == eCollisionBetweenTrac 286 if (fITStepStatus == eCollisionBetweenTracks)
>> 287 // if(fInteractionStep == false)
192 { 288 {
193 G4VITReactionProcess* pReactionProcess << 289 // TODO
194 fReactionInfo = pReactionProcess->Find << 290 FindReaction(fReactionSet,
195 currentTimeStep, << 291 currentTimeStep,
196 fGlobalTime, << 292 previousTimeStep,
197 reachedUserTimeLimit); << 293 reachedUserTimeLimit);
198 <<
199 // TODO 294 // TODO
200 // A ne faire uniquement si le temps c 295 // A ne faire uniquement si le temps choisis est celui calculé par le time stepper
201 // Sinon utiliser quelque chose comme 296 // Sinon utiliser quelque chose comme : fModelProcessor->FindReaction(&fMainList);
202 297
203 for (auto& pChanges : fReactionInfo) 298 for (auto& pChanges : fReactionInfo)
204 { 299 {
205 auto pTrackA = const_cast<G4Track* 300 auto pTrackA = const_cast<G4Track*>(pChanges->GetTrackA());
206 auto pTrackB = const_cast<G4Track* 301 auto pTrackB = const_cast<G4Track*>(pChanges->GetTrackB());
207 302
208 if (pTrackA == nullptr 303 if (pTrackA == nullptr
209 || pTrackB == nullptr 304 || pTrackB == nullptr
210 || pTrackA->GetTrackStatus() = 305 || pTrackA->GetTrackStatus() == fStopAndKill
211 || pTrackB->GetTrackStatus() = 306 || pTrackB->GetTrackStatus() == fStopAndKill)
212 { 307 {
213 continue; 308 continue;
214 } 309 }
215 310
216 G4int nbSecondaries = pChanges->Ge 311 G4int nbSecondaries = pChanges->GetNumberOfSecondaries();
217 const std::vector<G4Track*>* produ 312 const std::vector<G4Track*>* productsVector = pChanges->GetfSecondary();
218 313
219 if (fpUserTimeStepAction != nullpt << 314 if (fpUserTimeStepAction)
220 { 315 {
221 fpUserTimeStepAction->UserReac 316 fpUserTimeStepAction->UserReactionAction(*pTrackA,
222 317 *pTrackB,
223 318 productsVector);
224 } 319 }
225 320
226 #ifdef G4VERBOSE 321 #ifdef G4VERBOSE
227 if (fVerbose != 0) << 322 if (fVerbose)
228 { 323 {
229 G4cout << "At time : " << std: 324 G4cout << "At time : " << std::setw(7) << G4BestUnit(fGlobalTime, "Time")
230 << " Reaction : " << Ge 325 << " Reaction : " << GetIT(pTrackA)->GetName() << " ("
231 << pTrackA->GetTrackID( 326 << pTrackA->GetTrackID() << ") + " << GetIT(pTrackB)->GetName() << " ("
232 << pTrackB->GetTrackID( 327 << pTrackB->GetTrackID() << ") -> ";
233 } 328 }
234 #endif 329 #endif
235 330
236 if (nbSecondaries > 0) 331 if (nbSecondaries > 0)
237 { 332 {
238 for (int i = 0; i < nbSecondar 333 for (int i = 0; i < nbSecondaries; ++i)
239 { 334 {
240 #ifdef G4VERBOSE 335 #ifdef G4VERBOSE
241 if ((fVerbose != 0) && i ! << 336 if (fVerbose && i != 0)
242 { 337 {
243 G4cout << " + "; 338 G4cout << " + ";
244 } 339 }
245 #endif 340 #endif
246 341
247 G4Track* secondary = (*pro 342 G4Track* secondary = (*productsVector)[i]; //changes->GetSecondary(i);
248 // fpTrackContainer->_PushT << 343 fpTrackContainer->_PushTrack(secondary);
249 GetIT(secondary)->SetParen 344 GetIT(secondary)->SetParentID(pTrackA->GetTrackID(),
250 345 pTrackB->GetTrackID());
251 346
252 if (secondary->GetGlobalTi 347 if (secondary->GetGlobalTime() - fGlobalTime > fTimeTolerance)
253 { 348 {
254 G4ExceptionDescription 349 G4ExceptionDescription exceptionDescription;
255 exceptionDescription < 350 exceptionDescription << "The time of the secondary should not be bigger than the"
256 351 " current global time."
257 < 352 << " This may cause synchronization problem. If the process you"
258 353 " are using required "
259 < << 354 << "such feature please contact the developpers." << G4endl
260 < 355 << "The global time in the step manager : "
261 < 356 << G4BestUnit(fGlobalTime, "Time")
262 < 357 << G4endl
263 < 358 << "The global time of the track : "
264 < 359 << G4BestUnit(secondary->GetGlobalTime(), "Time")
265 < 360 << G4endl;
266 361
267 G4Exception("G4Schedul 362 G4Exception("G4Scheduler::ComputeInteractionBetweenTracks",
268 "ITSchedul 363 "ITScheduler010",
269 FatalError 364 FatalErrorInArgument,
270 exceptionD 365 exceptionDescription);
271 } 366 }
272 367
273 #ifdef G4VERBOSE 368 #ifdef G4VERBOSE
274 if (fVerbose != 0) << 369 if (fVerbose)
275 { 370 {
276 G4cout << GetIT(second 371 G4cout << GetIT(secondary)->GetName() << " ("
277 << secondary->G 372 << secondary->GetTrackID() << ")";
278 } 373 }
279 #endif 374 #endif
280 } 375 }
281 } 376 }
282 else 377 else
283 { 378 {
284 #ifdef G4VERBOSE 379 #ifdef G4VERBOSE
285 if (fVerbose != 0) << 380 if (fVerbose)
286 { 381 {
287 G4cout << "No product"; 382 G4cout << "No product";
288 } 383 }
289 #endif 384 #endif
290 } 385 }
291 #ifdef G4VERBOSE 386 #ifdef G4VERBOSE
292 if (fVerbose != 0) << 387 if (fVerbose)
293 { 388 {
294 G4cout << G4endl; 389 G4cout << G4endl;
295 } 390 }
296 #endif 391 #endif
297 if (pTrackA->GetTrackID() == 0 || 392 if (pTrackA->GetTrackID() == 0 || pTrackB->GetTrackID() == 0)
298 { 393 {
299 G4Track* pTrack = nullptr; 394 G4Track* pTrack = nullptr;
300 if (pTrackA->GetTrackID() == 0 395 if (pTrackA->GetTrackID() == 0)
301 { 396 {
302 pTrack = pTrackA; 397 pTrack = pTrackA;
303 } 398 }
304 else 399 else
305 { 400 {
306 pTrack = pTrackB; 401 pTrack = pTrackB;
307 } 402 }
308 403
309 G4ExceptionDescription excepti 404 G4ExceptionDescription exceptionDescription;
310 exceptionDescription 405 exceptionDescription
311 << "The problem was found 406 << "The problem was found for the reaction between tracks :"
312 << pTrackA->GetParticleDef 407 << pTrackA->GetParticleDefinition()->GetParticleName() << " ("
313 << pTrackA->GetTrackID() < 408 << pTrackA->GetTrackID() << ") & "
314 << pTrackB->GetParticleDef 409 << pTrackB->GetParticleDefinition()->GetParticleName() << " ("
315 << pTrackB->GetTrackID() < 410 << pTrackB->GetTrackID() << "). \n";
316 411
317 if (pTrack->GetStep() == nullp 412 if (pTrack->GetStep() == nullptr)
318 { 413 {
319 exceptionDescription << "A 414 exceptionDescription << "Also no step was found"
320 << " 415 << " ie track->GetStep() == 0 \n";
321 } 416 }
322 417
323 exceptionDescription << "Paren 418 exceptionDescription << "Parent ID of trackA : "
324 << pTrack 419 << pTrackA->GetParentID() << "\n";
325 exceptionDescription << "Paren 420 exceptionDescription << "Parent ID of trackB : "
326 << pTrack 421 << pTrackB->GetParentID() << "\n";
327 422
328 exceptionDescription 423 exceptionDescription
329 << "The ID of one of the r 424 << "The ID of one of the reaction track was not setup.";
330 G4Exception("G4Scheduler::Comp 425 G4Exception("G4Scheduler::ComputeInteractionBetweenTracks",
331 "ITScheduler011", 426 "ITScheduler011",
332 FatalErrorInArgume 427 FatalErrorInArgument,
333 exceptionDescripti 428 exceptionDescription);
334 } 429 }
335 430
336 if (pChanges->WereParentsKilled()) 431 if (pChanges->WereParentsKilled())
337 { 432 {
338 pTrackA->SetTrackStatus(fStopA 433 pTrackA->SetTrackStatus(fStopAndKill);
339 pTrackB->SetTrackStatus(fStopA 434 pTrackB->SetTrackStatus(fStopAndKill);
340 435
341 fpTrackingManager->EndTracking 436 fpTrackingManager->EndTracking(pTrackA);
342 fpTrackingManager->EndTracking 437 fpTrackingManager->EndTracking(pTrackB);
343 } 438 }
344 439
345 pChanges.reset(nullptr); 440 pChanges.reset(nullptr);
346 } 441 }
347 442
348 fReactionInfo.clear(); 443 fReactionInfo.clear();
349 } 444 }
350 445
351 // fReactionSet->CleanAllReaction(); << 446 fReactionSet->CleanAllReaction();
352 447
353 fpTrackContainer->MergeSecondariesWithMain 448 fpTrackContainer->MergeSecondariesWithMainList();
354 fpTrackContainer->KillTracks(); 449 fpTrackContainer->KillTracks();
>> 450 }
>> 451
>> 452 //______________________________________________________________________________
>> 453 void G4ITModelProcessor::FindReaction(G4ITReactionSet* pReactionSet,
>> 454 const double currentStepTime,
>> 455 const double /*previousStepTime*/,
>> 456 const bool reachedUserStepTimeLimit)
>> 457 {
>> 458 if (pReactionSet == nullptr || fActiveModels.empty())
>> 459 {
>> 460 return;
>> 461 }
>> 462
>> 463 G4ITReactionPerTrackMap& reactionPerTrackMap = pReactionSet->GetReactionMap();
>> 464 G4VITReactionProcess* pReactionProcess = fpActiveModelWithMinTimeStep->GetReactionProcess();
>> 465
>> 466 for (auto tracks_i = reactionPerTrackMap.begin();
>> 467 tracks_i != reactionPerTrackMap.end();
>> 468 tracks_i = reactionPerTrackMap.begin())
>> 469 {
>> 470 G4Track* pTrackA = tracks_i->first;
>> 471 if (pTrackA->GetTrackStatus() == fStopAndKill)
>> 472 {
>> 473 continue;
>> 474 }
>> 475
>> 476 G4ITReactionPerTrackPtr reactionPerTrack = tracks_i->second;
>> 477 G4ITReactionList& reactionList = reactionPerTrack->GetReactionList();
>> 478
>> 479 assert(reactionList.begin() != reactionList.end());
>> 480
>> 481 for (auto it = reactionList.begin(); it != reactionList.end(); it = reactionList.begin())
>> 482 {
>> 483 G4ITReactionPtr reaction(*it);
>> 484 G4Track* pTrackB = reaction->GetReactant(pTrackA);
>> 485 if (pTrackB->GetTrackStatus() == fStopAndKill)
>> 486 {
>> 487 continue;
>> 488 }
>> 489
>> 490 if (pTrackB == pTrackA)
>> 491 {
>> 492 G4ExceptionDescription exceptionDescription;
>> 493 exceptionDescription
>> 494 << "The IT reaction process sent back a reaction between trackA and trackB. ";
>> 495 exceptionDescription << "The problem is trackA == trackB";
>> 496 G4Exception("G4ITModelProcessor::FindReaction",
>> 497 "ITModelProcessor005",
>> 498 FatalErrorInArgument,
>> 499 exceptionDescription);
>> 500 }
>> 501
>> 502 pReactionSet->SelectThisReaction(reaction);
>> 503
>> 504 if (pReactionProcess && pReactionProcess->TestReactibility(*pTrackA,
>> 505 *pTrackB,
>> 506 currentStepTime,
>> 507 reachedUserStepTimeLimit))
>> 508 {
>> 509 auto pReactionChange = pReactionProcess->MakeReaction(*pTrackA, *pTrackB);
>> 510
>> 511 if (pReactionChange)
>> 512 {
>> 513 fReactionInfo.push_back(std::move(pReactionChange));
>> 514 break;
>> 515 }
>> 516 }
>> 517 }
>> 518 }
>> 519
>> 520 //assert(G4ITReaction::gAll->empty() == true);
355 } 521 }
356 522
357 void G4ITModelProcessor::SetTrack(const G4Trac 523 void G4ITModelProcessor::SetTrack(const G4Track* track)
358 { 524 {
359 fpTrack = track; 525 fpTrack = track;
360 } 526 }
361 527
362 void G4ITModelProcessor::SetModelHandler(G4ITM 528 void G4ITModelProcessor::SetModelHandler(G4ITModelHandler* pModelHandler)
363 { 529 {
364 if (fInitialized) 530 if (fInitialized)
365 { 531 {
366 G4ExceptionDescription exceptionDescri 532 G4ExceptionDescription exceptionDescription;
367 exceptionDescription 533 exceptionDescription
368 << "You are trying to set a new mo 534 << "You are trying to set a new model while the model processor has alreaday be initialized";
369 G4Exception("G4ITModelProcessor::SetMo 535 G4Exception("G4ITModelProcessor::SetModelHandler", "ITModelProcessor001",
370 FatalErrorInArgument, exce 536 FatalErrorInArgument, exceptionDescription);
371 } 537 }
372 fpModelHandler = pModelHandler; 538 fpModelHandler = pModelHandler;
373 } 539 }
374 540
375 void G4ITModelProcessor::CleanProcessor() 541 void G4ITModelProcessor::CleanProcessor()
376 { 542 {
377 fpTrack = nullptr; 543 fpTrack = nullptr;
378 } 544 }
379 545
380 bool G4ITModelProcessor::GetComputeTimeStep() 546 bool G4ITModelProcessor::GetComputeTimeStep() const
381 { 547 {
382 return fComputeTimeStep; 548 return fComputeTimeStep;
383 } 549 }
384 550
385 const G4Track* G4ITModelProcessor::GetTrack() 551 const G4Track* G4ITModelProcessor::GetTrack() const
386 { 552 {
387 return fpTrack; 553 return fpTrack;
388 } 554 }
389 555
390 void G4ITModelProcessor::SetTrackingManager(G4 556 void G4ITModelProcessor::SetTrackingManager(G4ITTrackingManager* pTrackingManager)
391 { 557 {
392 fpTrackingManager = pTrackingManager; 558 fpTrackingManager = pTrackingManager;
393 } 559 }
394 560
395 561