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