| 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 "G4VITTimeStepper.hh"
37 #include "G4VITReactionProcess.hh" 37 #include "G4VITReactionProcess.hh"
38 #include "G4ITReaction.hh" << 38
39 #include "G4ITTrackHolder.hh" << 39 std::map<const G4Track*, G4bool> G4ITModelProcessor::fHasReacted ;
40 #include "G4ITTrackingManager.hh" <<
41 #include "G4VITStepModel.hh" <<
42 #include "G4UserTimeStepAction.hh" <<
43 #include "G4UnitsTable.hh" <<
44 #include "G4Scheduler.hh" <<
45 #include "G4SystemOfUnits.hh" <<
46 #include <vector> <<
47 <<
48 //#define DEBUG_MEM <<
49 <<
50 #ifdef DEBUG_MEM <<
51 #include "G4MemStat.hh" <<
52 using namespace G4MemStat; <<
53 #endif <<
54 40
55 G4ITModelProcessor::G4ITModelProcessor() 41 G4ITModelProcessor::G4ITModelProcessor()
56 { 42 {
57 fpTrack = nullptr; << 43 //ctor
>> 44 fTrack = 0;
>> 45 fpModelHandler = 0;
>> 46 fpModel = 0;
58 fInitialized = false; 47 fInitialized = false;
59 fUserMinTimeStep = -1.; << 48 fpModelManager = 0;
60 fTSTimeStep = DBL_MAX; << 49 fCurrentModel.assign(G4ITType::size(), std::vector<G4VITModel*>());
61 fpTrackingManager = nullptr; << 50
62 fReactionSet = nullptr; << 51 for(int i = 0 ; i < (int) G4ITType::size() ; i++)
63 fpTrackContainer = nullptr; << 52 {
64 fpModelHandler = nullptr; << 53 fCurrentModel[i].assign(G4ITType::size(),0);
65 fpActiveModelWithMinTimeStep = nullptr; << 54 }
66 fComputeTimeStep = false; <<
67 fComputeReaction = false; <<
68 } 55 }
69 56
70 G4ITModelProcessor::~G4ITModelProcessor() = de << 57 G4ITModelProcessor::~G4ITModelProcessor()
>> 58 {
>> 59 //dtor
>> 60 if(fpModelHandler) delete fpModelHandler;
>> 61 fCurrentModel.clear();
>> 62 fReactionInfo.clear();
>> 63 }
71 64
72 void G4ITModelProcessor::RegisterModel(double << 65 // Should not be used
>> 66 G4ITModelProcessor::G4ITModelProcessor(const G4ITModelProcessor& /*other*/)
73 { 67 {
74 fpModelHandler->RegisterModel(model, time) << 68 //copy ctorr
>> 69 fTrack = 0;
>> 70 fpModelHandler = 0;
>> 71 fpModel = 0;
>> 72 fInitialized = false;
>> 73 fpModelManager = 0;
75 } 74 }
76 75
>> 76 // Should not be used
>> 77 G4ITModelProcessor& G4ITModelProcessor::operator=(const G4ITModelProcessor& rhs)
>> 78 {
>> 79 if (this == &rhs) return *this; // handle self assignment
>> 80 //assignment operator
>> 81 return *this;
>> 82 }
>> 83 //______________________________________________________________________________
77 void G4ITModelProcessor::Initialize() 84 void G4ITModelProcessor::Initialize()
78 { 85 {
79 fpModelHandler->Initialize(); 86 fpModelHandler->Initialize();
80 fReactionSet = G4ITReactionSet::Instance() <<
81 fpTrackContainer = G4ITTrackHolder::Instan <<
82 fInitialized = true; 87 fInitialized = true;
83 fComputeTimeStep = false; << 88 }
84 fComputeReaction = false; << 89
85 if (fpModelHandler->GetTimeStepComputerFla << 90 //______________________________________________________________________________
>> 91 void G4ITModelProcessor::InitializeStepper(const G4double& currentGlobalTime,
>> 92 const G4double& userMinTime)
>> 93 {
>> 94 // G4cout << "G4ITModelProcessor::InitializeStepper" << G4endl;
>> 95 if(fpModelHandler==0)
86 { 96 {
87 fComputeTimeStep = true; << 97 G4ExceptionDescription exceptionDescription ;
>> 98 exceptionDescription << "No G4ITModelHandler was passed to the modelProcessor.";
>> 99 G4Exception("G4ITModelProcessor::InitializeStepper","ITModelProcessor002",
>> 100 FatalErrorInArgument,exceptionDescription);
88 } 101 }
89 if (fpModelHandler->GetReactionProcessFlag << 102 const std::vector<std::vector<G4ITModelManager*> >* modelManager = fpModelHandler
>> 103 ->GetAllModelManager();
>> 104
>> 105 if(modelManager==0)
90 { 106 {
91 fComputeReaction = true; << 107 G4ExceptionDescription exceptionDescription ;
>> 108 exceptionDescription << "No G4ITModelManager was register to G4ITModelHandler.";
>> 109 G4Exception("G4ITModelProcessor::InitializeStepper","ITModelProcessor003",
>> 110 FatalErrorInArgument,exceptionDescription);
92 } 111 }
93 } <<
94 112
95 G4double G4ITModelProcessor::CalculateMinTimeS << 113 int nbModels1 = modelManager->size() ;
96 <<
97 { <<
98 114
99 #if defined (DEBUG_MEM) && defined (DEBUG_MEM_ << 115 G4VITTimeStepper::SetTimes(currentGlobalTime, userMinTime) ;
100 MemStat mem_first, mem_second, mem_diff; <<
101 mem_first = MemoryUsage(); <<
102 #endif <<
103 <<
104 fpActiveModelWithMinTimeStep = nullptr; <<
105 fTSTimeStep = DBL_MAX; <<
106 <<
107 InitializeStepper(currentGlobalTime, defin <<
108 <<
109 #if defined (DEBUG_MEM) && defined (DEBUG_MEM_ <<
110 mem_second = MemoryUsage(); <<
111 mem_diff = mem_second-mem_first; <<
112 G4cout << "\t || MEM || G4Scheduler::Calcu <<
113 "computing fpModelProcessor -> InitializeS <<
114 << mem_diff <<
115 << G4endl; <<
116 #endif <<
117 116
118 for (auto& pStepModel : fActiveModels) << 117 // TODO !!!
>> 118 // if( nbModels1 != 1 || (nbModels1 == 1 && !fpModelManager) )
119 { 119 {
120 fTSTimeStep = << 120 int nbModels2 = -1;
121 pStepModel->GetTimeStepper()->Calc << 121 G4VITModel* model = 0;
122 currentGlobalTime, << 122 G4ITModelManager* modman = 0;
123 definedMinTimeStep); << 123
124 << 124 for(int i = 0 ; i < nbModels1 ; i++)
125 fpActiveModelWithMinTimeStep = pStepMo << 125 {
126 << 126 nbModels2 = (*modelManager)[i].size();
127 if(fTSTimeStep == -1){ << 127
128 fpActiveModelWithMinTimeStep->GetR << 128 for(int j = 0 ; j <= i ; j++)
129 if(fReactionSet->Empty()) return D << 129 {
130 const auto& fReactionSetInTime = f << 130 modman = (*modelManager)[i][j];
131 fTSTimeStep = fReactionSetInTime.b << 131
>> 132 if(modman == 0) continue ;
>> 133
>> 134 model = modman -> GetModel(currentGlobalTime);
>> 135 G4VITTimeStepper* stepper = model->GetTimeStepper() ;
>> 136
>> 137 stepper -> PrepareForAllProcessors() ;
>> 138 stepper -> Prepare() ;
>> 139 fCurrentModel[i][j] = model;
>> 140 }
>> 141 }
>> 142
>> 143 if(nbModels1 == 1 && nbModels2 ==1)
>> 144 {
>> 145 fpModelManager = modman;
>> 146 fpModel = model;
132 } 147 }
>> 148 else fpModel = 0;
>> 149 }
>> 150 }
>> 151
>> 152 //______________________________________________________________________________
>> 153 void G4ITModelProcessor::CalculateTimeStep(const G4Track* track, const G4double userMinTimeStep)
>> 154 {
>> 155 // G4cout << "G4ITModelProcessor::CalculateStep" << G4endl;
>> 156 CleanProcessor();
>> 157 if(track == 0)
>> 158 {
>> 159 G4ExceptionDescription exceptionDescription ;
>> 160 exceptionDescription << "No track was passed to the method (track == 0).";
>> 161 G4Exception("G4ITModelProcessor::CalculateStep","ITModelProcessor004",
>> 162 FatalErrorInArgument,exceptionDescription);
133 } 163 }
>> 164 SetTrack(track);
>> 165 fUserMinTimeStep = userMinTimeStep ;
134 166
135 #if defined (DEBUG_MEM) && defined (DEBUG_MEM_ << 167 DoCalculateStep();
136 mem_second = MemoryUsage(); <<
137 mem_diff = mem_second-mem_first; <<
138 G4cout << "\t || MEM || G4Scheduler::Calcu <<
139 "After looping on tracks, diff is : " << m <<
140 #endif <<
141 return fTSTimeStep; <<
142 } 168 }
143 169
144 //____________________________________________ 170 //______________________________________________________________________________
>> 171 void G4ITModelProcessor::DoCalculateStep()
>> 172 {
>> 173 if(fpModel) // ie only one model has been declared and will be used
>> 174 {
>> 175 fpModel -> GetTimeStepper()->CalculateStep(*fTrack, fUserMinTimeStep);
>> 176 }
>> 177 else // ie many models have been declared and will be used
>> 178 {
>> 179 std::vector<G4VITModel*>& model = fCurrentModel[GetIT(fTrack)->GetITType()];
145 180
146 void G4ITModelProcessor::InitializeStepper(G4d << 181 for(int i =0 ; i < (int) model.size() ; i++)
147 G4d << 182 {
>> 183 if(model[i] == 0) continue;
>> 184 model[i]->GetTimeStepper()->CalculateStep(*fTrack, fUserMinTimeStep);
>> 185 }
>> 186 }
>> 187 }
>> 188
>> 189 //______________________________________________________________________________
>> 190 void G4ITModelProcessor::FindReaction(std::map<G4Track*, G4TrackVectorHandle>* tracks,
>> 191 const double currentStepTime,
>> 192 const double previousStepTime,
>> 193 const bool reachedUserStepTimeLimit)
148 { 194 {
149 G4VITTimeStepComputer::SetTimes(currentGlo << 195 // DEBUG
>> 196 // G4cout << "G4ITReactionManager::FindReaction" << G4endl;
>> 197 if(tracks == 0) return ;
150 198
151 #if defined (DEBUG_MEM) << 199 if(fpModelHandler->GetAllModelManager()->empty()) return ;
152 MemStat mem_first, mem_second, mem_diff; <<
153 mem_first = MemoryUsage(); <<
154 #endif <<
155 200
156 fActiveModels = fpModelHandler->GetActiveM << 201 std::map<G4Track*, G4TrackVectorHandle>::iterator tracks_i = tracks->begin();;
157 202
158 for (auto& pModel : fActiveModels) << 203 for(tracks_i = tracks->begin() ; tracks_i != tracks-> end() ; tracks_i ++)
159 { 204 {
160 pModel->PrepareNewTimeStep(); << 205 /// Get track A
161 } << 206 G4Track* trackA = tracks_i->first;
162 207
163 #if defined (DEBUG_MEM) << 208 if(trackA == 0) continue;
164 mem_second = MemoryUsage(); <<
165 mem_diff = mem_second-mem_first; <<
166 G4cout << "\t || MEM || G4ITModelP <<
167 #endif <<
168 209
169 } << 210 std::map<const G4Track*, G4bool>::iterator it_hasReacted = fHasReacted.find(trackA);
>> 211 if(it_hasReacted != fHasReacted.end()) continue;
>> 212 if(trackA->GetTrackStatus() == fStopAndKill) continue;
170 213
171 //____________________________________________ << 214 G4IT* ITA = GetIT(trackA);
>> 215 G4ITType ITypeA = ITA -> GetITType();
172 216
173 void G4ITModelProcessor::ComputeTrackReaction( << 217 const std::vector<G4VITModel*> model = fCurrentModel[ITypeA];
174 <<
175 <<
176 <<
177 <<
178 <<
179 <<
180 <<
181 #ifdef G4VERBOSE <<
182 fVerbose <<
183 #endif <<
184 ) <<
185 { <<
186 if (fReactionSet->Empty()) <<
187 { <<
188 return; <<
189 } <<
190 218
191 if (fITStepStatus == eCollisionBetweenTrac << 219 G4TrackVectorHandle& trackB_vector = tracks_i->second ;
192 { << 220 std::vector<G4Track*>::iterator trackB_i = trackB_vector->begin();
193 G4VITReactionProcess* pReactionProcess <<
194 fReactionInfo = pReactionProcess->Find <<
195 currentTimeStep, <<
196 fGlobalTime, <<
197 reachedUserTimeLimit); <<
198 <<
199 // TODO <<
200 // A ne faire uniquement si le temps c <<
201 // Sinon utiliser quelque chose comme <<
202 221
203 for (auto& pChanges : fReactionInfo) << 222 G4Track* trackB = 0 ;
>> 223 G4ITType ITypeB(-1);
>> 224 G4VITReactionProcess* process = 0;
>> 225 G4ITReactionChange* changes = 0;
>> 226
>> 227 for(; trackB_i != trackB_vector->end() ; trackB_i++)
204 { 228 {
205 auto pTrackA = const_cast<G4Track* << 229 trackB = *trackB_i;
206 auto pTrackB = const_cast<G4Track* <<
207 230
208 if (pTrackA == nullptr << 231 if(trackB == 0) continue;
209 || pTrackB == nullptr << 232 it_hasReacted = fHasReacted.find(trackB);
210 || pTrackA->GetTrackStatus() = << 233 if(it_hasReacted != fHasReacted.end()) continue;
211 || pTrackB->GetTrackStatus() = << 234 if(trackB->GetTrackStatus() == fStopAndKill) continue;
212 { <<
213 continue; <<
214 } <<
215 235
216 G4int nbSecondaries = pChanges->Ge << 236 // DEBUG
217 const std::vector<G4Track*>* produ << 237 // G4cout << "Couple : " << trackA->GetParticleDefinition->GetParticleName() << " ("
>> 238 // << trackA->GetTrackID() << ") "
>> 239 // << trackB->GetParticleDefinition->GetParticleName() << " ("
>> 240 // << trackB->GetTrackID() << ")"
>> 241 // << G4endl;
218 242
219 if (fpUserTimeStepAction != nullpt << 243 if(trackB == trackA)
220 { 244 {
221 fpUserTimeStepAction->UserReac << 245 G4ExceptionDescription exceptionDescription ;
222 << 246 exceptionDescription << "The IT reaction process sent back a reaction between trackA and trackB. ";
223 << 247 exceptionDescription << "The problem is trackA == trackB";
>> 248 G4Exception("G4ITModelProcessor::FindReaction","ITModelProcessor005",
>> 249 FatalErrorInArgument,exceptionDescription);
224 } 250 }
225 251
226 #ifdef G4VERBOSE << 252 G4IT* ITB = GetIT(trackB);
227 if (fVerbose != 0) << 253 G4ITType ITypeBtmp = ITB -> GetITType();
228 { <<
229 G4cout << "At time : " << std: <<
230 << " Reaction : " << Ge <<
231 << pTrackA->GetTrackID( <<
232 << pTrackB->GetTrackID( <<
233 } <<
234 #endif <<
235 254
236 if (nbSecondaries > 0) << 255 if(ITypeB != ITypeBtmp)
237 { 256 {
238 for (int i = 0; i < nbSecondar << 257 ITypeB = ITypeBtmp ;
239 { << 258
240 #ifdef G4VERBOSE << 259 if(model[ITypeB])
241 if ((fVerbose != 0) && i ! << 260 process = model[ITypeB]->GetReactionProcess();
242 { <<
243 G4cout << " + "; <<
244 } <<
245 #endif <<
246 <<
247 G4Track* secondary = (*pro <<
248 // fpTrackContainer->_PushT <<
249 GetIT(secondary)->SetParen <<
250 <<
251 <<
252 if (secondary->GetGlobalTi <<
253 { <<
254 G4ExceptionDescription <<
255 exceptionDescription < <<
256 <<
257 < <<
258 <<
259 < <<
260 < <<
261 < <<
262 < <<
263 < <<
264 < <<
265 < <<
266 <<
267 G4Exception("G4Schedul <<
268 "ITSchedul <<
269 FatalError <<
270 exceptionD <<
271 } <<
272 <<
273 #ifdef G4VERBOSE <<
274 if (fVerbose != 0) <<
275 { <<
276 G4cout << GetIT(second <<
277 << secondary->G <<
278 } <<
279 #endif <<
280 } <<
281 } <<
282 else <<
283 { <<
284 #ifdef G4VERBOSE <<
285 if (fVerbose != 0) <<
286 { <<
287 G4cout << "No product"; <<
288 } <<
289 #endif <<
290 } <<
291 #ifdef G4VERBOSE <<
292 if (fVerbose != 0) <<
293 { <<
294 G4cout << G4endl; <<
295 } <<
296 #endif <<
297 if (pTrackA->GetTrackID() == 0 || <<
298 { <<
299 G4Track* pTrack = nullptr; <<
300 if (pTrackA->GetTrackID() == 0 <<
301 { <<
302 pTrack = pTrackA; <<
303 } <<
304 else <<
305 { <<
306 pTrack = pTrackB; <<
307 } <<
308 <<
309 G4ExceptionDescription excepti <<
310 exceptionDescription <<
311 << "The problem was found <<
312 << pTrackA->GetParticleDef <<
313 << pTrackA->GetTrackID() < <<
314 << pTrackB->GetParticleDef <<
315 << pTrackB->GetTrackID() < <<
316 <<
317 if (pTrack->GetStep() == nullp <<
318 { <<
319 exceptionDescription << "A <<
320 << " <<
321 } <<
322 <<
323 exceptionDescription << "Paren <<
324 << pTrack <<
325 exceptionDescription << "Paren <<
326 << pTrack <<
327 <<
328 exceptionDescription <<
329 << "The ID of one of the r <<
330 G4Exception("G4Scheduler::Comp <<
331 "ITScheduler011", <<
332 FatalErrorInArgume <<
333 exceptionDescripti <<
334 } 261 }
335 262
336 if (pChanges->WereParentsKilled()) << 263 if(process && process -> TestReactibility(*trackA, *trackB,
>> 264 currentStepTime, previousStepTime,
>> 265 reachedUserStepTimeLimit))
337 { 266 {
338 pTrackA->SetTrackStatus(fStopA << 267 changes = process->MakeReaction(*trackA, *trackB);
339 pTrackB->SetTrackStatus(fStopA <<
340 <<
341 fpTrackingManager->EndTracking <<
342 fpTrackingManager->EndTracking <<
343 } 268 }
344 269
345 pChanges.reset(nullptr); << 270 if(changes)
346 } << 271 {
347 << 272 fHasReacted[trackA] = true;
348 fReactionInfo.clear(); << 273 fHasReacted[trackB] = true;
349 } << 274 changes -> GetTrackA();
350 << 275 changes -> GetTrackB();
351 // fReactionSet->CleanAllReaction(); <<
352 276
353 fpTrackContainer->MergeSecondariesWithMain << 277 fReactionInfo.push_back(changes);
354 fpTrackContainer->KillTracks(); <<
355 } <<
356 278
357 void G4ITModelProcessor::SetTrack(const G4Trac << 279 process->ResetChanges();
358 { << 280 changes = 0;
359 fpTrack = track; <<
360 } <<
361 281
362 void G4ITModelProcessor::SetModelHandler(G4ITM << 282 break;
363 { << 283 }
364 if (fInitialized) << 284 }
365 { <<
366 G4ExceptionDescription exceptionDescri <<
367 exceptionDescription <<
368 << "You are trying to set a new mo <<
369 G4Exception("G4ITModelProcessor::SetMo <<
370 FatalErrorInArgument, exce <<
371 } 285 }
372 fpModelHandler = pModelHandler; <<
373 } <<
374 286
375 void G4ITModelProcessor::CleanProcessor() << 287 fHasReacted.clear();
376 { <<
377 fpTrack = nullptr; <<
378 } 288 }
379 <<
380 bool G4ITModelProcessor::GetComputeTimeStep() <<
381 { <<
382 return fComputeTimeStep; <<
383 } <<
384 <<
385 const G4Track* G4ITModelProcessor::GetTrack() <<
386 { <<
387 return fpTrack; <<
388 } <<
389 <<
390 void G4ITModelProcessor::SetTrackingManager(G4 <<
391 { <<
392 fpTrackingManager = pTrackingManager; <<
393 } <<
394 <<
395 289