| 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& pStepModel : fActiveModels)
119 { 117 {
120 fTSTimeStep = 118 fTSTimeStep =
121 pStepModel->GetTimeStepper()->Calc 119 pStepModel->GetTimeStepper()->CalculateMinTimeStep(
122 currentGlobalTime, 120 currentGlobalTime,
123 definedMinTimeStep); 121 definedMinTimeStep);
124 122
125 fpActiveModelWithMinTimeStep = pStepMo 123 fpActiveModelWithMinTimeStep = pStepModel;
126 124
127 if(fTSTimeStep == -1){ 125 if(fTSTimeStep == -1){
128 fpActiveModelWithMinTimeStep->GetR 126 fpActiveModelWithMinTimeStep->GetReactionProcess()->Initialize();
129 if(fReactionSet->Empty()) return D 127 if(fReactionSet->Empty()) return DBL_MAX;
130 const auto& fReactionSetInTime = f << 128 auto fReactionSetInTime = fReactionSet->GetReactionsPerTime();
131 fTSTimeStep = fReactionSetInTime.b 129 fTSTimeStep = fReactionSetInTime.begin()->get()->GetTime() - currentGlobalTime;
132 } 130 }
133 } 131 }
134 132
135 #if defined (DEBUG_MEM) && defined (DEBUG_MEM_ 133 #if defined (DEBUG_MEM) && defined (DEBUG_MEM_DETAILED_STEPPING)
136 mem_second = MemoryUsage(); 134 mem_second = MemoryUsage();
137 mem_diff = mem_second-mem_first; 135 mem_diff = mem_second-mem_first;
138 G4cout << "\t || MEM || G4Scheduler::Calcu 136 G4cout << "\t || MEM || G4Scheduler::CalculateMinTimeStep || "
139 "After looping on tracks, diff is : " << m 137 "After looping on tracks, diff is : " << mem_diff << G4endl;
140 #endif 138 #endif
141 return fTSTimeStep; 139 return fTSTimeStep;
142 } 140 }
143 141
144 //____________________________________________ 142 //______________________________________________________________________________
145 143
146 void G4ITModelProcessor::InitializeStepper(G4d 144 void G4ITModelProcessor::InitializeStepper(G4double currentGlobalTime,
147 G4d 145 G4double userMinTime)
148 { 146 {
149 G4VITTimeStepComputer::SetTimes(currentGlo 147 G4VITTimeStepComputer::SetTimes(currentGlobalTime, userMinTime);
150 148
151 #if defined (DEBUG_MEM) 149 #if defined (DEBUG_MEM)
152 MemStat mem_first, mem_second, mem_diff; 150 MemStat mem_first, mem_second, mem_diff;
153 mem_first = MemoryUsage(); 151 mem_first = MemoryUsage();
154 #endif 152 #endif
155 153
156 fActiveModels = fpModelHandler->GetActiveM 154 fActiveModels = fpModelHandler->GetActiveModels(currentGlobalTime);
157 155
158 for (auto& pModel : fActiveModels) 156 for (auto& pModel : fActiveModels)
159 { 157 {
160 pModel->PrepareNewTimeStep(); 158 pModel->PrepareNewTimeStep();
161 } 159 }
162 160
163 #if defined (DEBUG_MEM) 161 #if defined (DEBUG_MEM)
164 mem_second = MemoryUsage(); 162 mem_second = MemoryUsage();
165 mem_diff = mem_second-mem_first; 163 mem_diff = mem_second-mem_first;
166 G4cout << "\t || MEM || G4ITModelP 164 G4cout << "\t || MEM || G4ITModelProcessor::InitializeStepper || After computing stepper -> Prepare(), diff is : " << mem_diff << G4endl;
167 #endif 165 #endif
168 166
169 } 167 }
170 168
171 //____________________________________________ 169 //_________________________________________________________________________
172 170
173 void G4ITModelProcessor::ComputeTrackReaction( 171 void G4ITModelProcessor::ComputeTrackReaction(G4ITStepStatus fITStepStatus,
174 172 G4double fGlobalTime,
175 173 G4double currentTimeStep,
176 174 G4double /*previousTimeStep*/,
177 175 G4bool reachedUserTimeLimit,
178 176 G4double fTimeTolerance,
179 177 G4UserTimeStepAction* fpUserTimeStepAction,
180 178 G4int
181 #ifdef G4VERBOSE 179 #ifdef G4VERBOSE
182 fVerbose 180 fVerbose
183 #endif 181 #endif
184 ) 182 )
185 { 183 {
186 if (fReactionSet->Empty()) 184 if (fReactionSet->Empty())
187 { 185 {
188 return; 186 return;
189 } 187 }
190 188
191 if (fITStepStatus == eCollisionBetweenTrac 189 if (fITStepStatus == eCollisionBetweenTracks)
192 { 190 {
193 G4VITReactionProcess* pReactionProcess 191 G4VITReactionProcess* pReactionProcess = fpActiveModelWithMinTimeStep->GetReactionProcess();
194 fReactionInfo = pReactionProcess->Find 192 fReactionInfo = pReactionProcess->FindReaction(fReactionSet,
195 currentTimeStep, 193 currentTimeStep,
196 fGlobalTime, 194 fGlobalTime,
197 reachedUserTimeLimit); 195 reachedUserTimeLimit);
198 196
199 // TODO 197 // TODO
200 // A ne faire uniquement si le temps c 198 // A ne faire uniquement si le temps choisis est celui calculé par le time stepper
201 // Sinon utiliser quelque chose comme 199 // Sinon utiliser quelque chose comme : fModelProcessor->FindReaction(&fMainList);
202 200
203 for (auto& pChanges : fReactionInfo) 201 for (auto& pChanges : fReactionInfo)
204 { 202 {
205 auto pTrackA = const_cast<G4Track* 203 auto pTrackA = const_cast<G4Track*>(pChanges->GetTrackA());
206 auto pTrackB = const_cast<G4Track* 204 auto pTrackB = const_cast<G4Track*>(pChanges->GetTrackB());
207 205
208 if (pTrackA == nullptr 206 if (pTrackA == nullptr
209 || pTrackB == nullptr 207 || pTrackB == nullptr
210 || pTrackA->GetTrackStatus() = 208 || pTrackA->GetTrackStatus() == fStopAndKill
211 || pTrackB->GetTrackStatus() = 209 || pTrackB->GetTrackStatus() == fStopAndKill)
212 { 210 {
213 continue; 211 continue;
214 } 212 }
215 213
216 G4int nbSecondaries = pChanges->Ge 214 G4int nbSecondaries = pChanges->GetNumberOfSecondaries();
217 const std::vector<G4Track*>* produ 215 const std::vector<G4Track*>* productsVector = pChanges->GetfSecondary();
218 216
219 if (fpUserTimeStepAction != nullpt << 217 if (fpUserTimeStepAction)
220 { 218 {
221 fpUserTimeStepAction->UserReac 219 fpUserTimeStepAction->UserReactionAction(*pTrackA,
222 220 *pTrackB,
223 221 productsVector);
224 } 222 }
225 223
226 #ifdef G4VERBOSE 224 #ifdef G4VERBOSE
227 if (fVerbose != 0) << 225 if (fVerbose)
228 { 226 {
229 G4cout << "At time : " << std: 227 G4cout << "At time : " << std::setw(7) << G4BestUnit(fGlobalTime, "Time")
230 << " Reaction : " << Ge 228 << " Reaction : " << GetIT(pTrackA)->GetName() << " ("
231 << pTrackA->GetTrackID( 229 << pTrackA->GetTrackID() << ") + " << GetIT(pTrackB)->GetName() << " ("
232 << pTrackB->GetTrackID( 230 << pTrackB->GetTrackID() << ") -> ";
233 } 231 }
234 #endif 232 #endif
235 233
236 if (nbSecondaries > 0) 234 if (nbSecondaries > 0)
237 { 235 {
238 for (int i = 0; i < nbSecondar 236 for (int i = 0; i < nbSecondaries; ++i)
239 { 237 {
240 #ifdef G4VERBOSE 238 #ifdef G4VERBOSE
241 if ((fVerbose != 0) && i ! << 239 if (fVerbose && i != 0)
242 { 240 {
243 G4cout << " + "; 241 G4cout << " + ";
244 } 242 }
245 #endif 243 #endif
246 244
247 G4Track* secondary = (*pro 245 G4Track* secondary = (*productsVector)[i]; //changes->GetSecondary(i);
248 // fpTrackContainer->_PushT 246 // fpTrackContainer->_PushTrack(secondary);
249 GetIT(secondary)->SetParen 247 GetIT(secondary)->SetParentID(pTrackA->GetTrackID(),
250 248 pTrackB->GetTrackID());
251 249
252 if (secondary->GetGlobalTi 250 if (secondary->GetGlobalTime() - fGlobalTime > fTimeTolerance)
253 { 251 {
254 G4ExceptionDescription 252 G4ExceptionDescription exceptionDescription;
255 exceptionDescription < 253 exceptionDescription << "The time of the secondary should not be bigger than the"
256 254 " current global time."
257 < 255 << " This may cause synchronization problem. If the process you"
258 256 " are using required "
259 < 257 << "such feature please contact the developers." << G4endl
260 < 258 << "The global time in the step manager : "
261 < 259 << G4BestUnit(fGlobalTime, "Time")
262 < 260 << G4endl
263 < 261 << "The global time of the track : "
264 < 262 << G4BestUnit(secondary->GetGlobalTime(), "Time")
265 < 263 << G4endl;
266 264
267 G4Exception("G4Schedul 265 G4Exception("G4Scheduler::ComputeInteractionBetweenTracks",
268 "ITSchedul 266 "ITScheduler010",
269 FatalError 267 FatalErrorInArgument,
270 exceptionD 268 exceptionDescription);
271 } 269 }
272 270
273 #ifdef G4VERBOSE 271 #ifdef G4VERBOSE
274 if (fVerbose != 0) << 272 if (fVerbose)
275 { 273 {
276 G4cout << GetIT(second 274 G4cout << GetIT(secondary)->GetName() << " ("
277 << secondary->G 275 << secondary->GetTrackID() << ")";
278 } 276 }
279 #endif 277 #endif
280 } 278 }
281 } 279 }
282 else 280 else
283 { 281 {
284 #ifdef G4VERBOSE 282 #ifdef G4VERBOSE
285 if (fVerbose != 0) << 283 if (fVerbose)
286 { 284 {
287 G4cout << "No product"; 285 G4cout << "No product";
288 } 286 }
289 #endif 287 #endif
290 } 288 }
291 #ifdef G4VERBOSE 289 #ifdef G4VERBOSE
292 if (fVerbose != 0) << 290 if (fVerbose)
293 { 291 {
294 G4cout << G4endl; 292 G4cout << G4endl;
295 } 293 }
296 #endif 294 #endif
297 if (pTrackA->GetTrackID() == 0 || 295 if (pTrackA->GetTrackID() == 0 || pTrackB->GetTrackID() == 0)
298 { 296 {
299 G4Track* pTrack = nullptr; 297 G4Track* pTrack = nullptr;
300 if (pTrackA->GetTrackID() == 0 298 if (pTrackA->GetTrackID() == 0)
301 { 299 {
302 pTrack = pTrackA; 300 pTrack = pTrackA;
303 } 301 }
304 else 302 else
305 { 303 {
306 pTrack = pTrackB; 304 pTrack = pTrackB;
307 } 305 }
308 306
309 G4ExceptionDescription excepti 307 G4ExceptionDescription exceptionDescription;
310 exceptionDescription 308 exceptionDescription
311 << "The problem was found 309 << "The problem was found for the reaction between tracks :"
312 << pTrackA->GetParticleDef 310 << pTrackA->GetParticleDefinition()->GetParticleName() << " ("
313 << pTrackA->GetTrackID() < 311 << pTrackA->GetTrackID() << ") & "
314 << pTrackB->GetParticleDef 312 << pTrackB->GetParticleDefinition()->GetParticleName() << " ("
315 << pTrackB->GetTrackID() < 313 << pTrackB->GetTrackID() << "). \n";
316 314
317 if (pTrack->GetStep() == nullp 315 if (pTrack->GetStep() == nullptr)
318 { 316 {
319 exceptionDescription << "A 317 exceptionDescription << "Also no step was found"
320 << " 318 << " ie track->GetStep() == 0 \n";
321 } 319 }
322 320
323 exceptionDescription << "Paren 321 exceptionDescription << "Parent ID of trackA : "
324 << pTrack 322 << pTrackA->GetParentID() << "\n";
325 exceptionDescription << "Paren 323 exceptionDescription << "Parent ID of trackB : "
326 << pTrack 324 << pTrackB->GetParentID() << "\n";
327 325
328 exceptionDescription 326 exceptionDescription
329 << "The ID of one of the r 327 << "The ID of one of the reaction track was not setup.";
330 G4Exception("G4Scheduler::Comp 328 G4Exception("G4Scheduler::ComputeInteractionBetweenTracks",
331 "ITScheduler011", 329 "ITScheduler011",
332 FatalErrorInArgume 330 FatalErrorInArgument,
333 exceptionDescripti 331 exceptionDescription);
334 } 332 }
335 333
336 if (pChanges->WereParentsKilled()) 334 if (pChanges->WereParentsKilled())
337 { 335 {
338 pTrackA->SetTrackStatus(fStopA 336 pTrackA->SetTrackStatus(fStopAndKill);
339 pTrackB->SetTrackStatus(fStopA 337 pTrackB->SetTrackStatus(fStopAndKill);
340 338
341 fpTrackingManager->EndTracking 339 fpTrackingManager->EndTracking(pTrackA);
342 fpTrackingManager->EndTracking 340 fpTrackingManager->EndTracking(pTrackB);
343 } 341 }
344 342
345 pChanges.reset(nullptr); 343 pChanges.reset(nullptr);
346 } 344 }
347 345
348 fReactionInfo.clear(); 346 fReactionInfo.clear();
349 } 347 }
350 348
351 // fReactionSet->CleanAllReaction(); 349 // fReactionSet->CleanAllReaction();
352 350
353 fpTrackContainer->MergeSecondariesWithMain 351 fpTrackContainer->MergeSecondariesWithMainList();
354 fpTrackContainer->KillTracks(); 352 fpTrackContainer->KillTracks();
355 } 353 }
356 354
357 void G4ITModelProcessor::SetTrack(const G4Trac 355 void G4ITModelProcessor::SetTrack(const G4Track* track)
358 { 356 {
359 fpTrack = track; 357 fpTrack = track;
360 } 358 }
361 359
362 void G4ITModelProcessor::SetModelHandler(G4ITM 360 void G4ITModelProcessor::SetModelHandler(G4ITModelHandler* pModelHandler)
363 { 361 {
364 if (fInitialized) 362 if (fInitialized)
365 { 363 {
366 G4ExceptionDescription exceptionDescri 364 G4ExceptionDescription exceptionDescription;
367 exceptionDescription 365 exceptionDescription
368 << "You are trying to set a new mo 366 << "You are trying to set a new model while the model processor has alreaday be initialized";
369 G4Exception("G4ITModelProcessor::SetMo 367 G4Exception("G4ITModelProcessor::SetModelHandler", "ITModelProcessor001",
370 FatalErrorInArgument, exce 368 FatalErrorInArgument, exceptionDescription);
371 } 369 }
372 fpModelHandler = pModelHandler; 370 fpModelHandler = pModelHandler;
373 } 371 }
374 372
375 void G4ITModelProcessor::CleanProcessor() 373 void G4ITModelProcessor::CleanProcessor()
376 { 374 {
377 fpTrack = nullptr; 375 fpTrack = nullptr;
378 } 376 }
379 377
380 bool G4ITModelProcessor::GetComputeTimeStep() 378 bool G4ITModelProcessor::GetComputeTimeStep() const
381 { 379 {
382 return fComputeTimeStep; 380 return fComputeTimeStep;
383 } 381 }
384 382
385 const G4Track* G4ITModelProcessor::GetTrack() 383 const G4Track* G4ITModelProcessor::GetTrack() const
386 { 384 {
387 return fpTrack; 385 return fpTrack;
388 } 386 }
389 387
390 void G4ITModelProcessor::SetTrackingManager(G4 388 void G4ITModelProcessor::SetTrackingManager(G4ITTrackingManager* pTrackingManager)
391 { 389 {
392 fpTrackingManager = pTrackingManager; 390 fpTrackingManager = pTrackingManager;
393 } 391 }
394 392
395 393