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