| 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 "G4ITTimeStepper.hh"
38 #include "G4ITReaction.hh" 39 #include "G4ITReaction.hh"
39 #include "G4ITTrackHolder.hh" 40 #include "G4ITTrackHolder.hh"
40 #include "G4ITTrackingManager.hh" 41 #include "G4ITTrackingManager.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
>> 46 using namespace std;
>> 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 = 0;
58 fInitialized = false; << 58 fpModel = 0;
59 fUserMinTimeStep = -1.; << 59 fInitialized = false;
60 fTSTimeStep = DBL_MAX; << 60 fpModelManager = 0;
61 fpTrackingManager = nullptr; << 61 fCurrentModel.assign(G4ITType::size(), std::vector<G4VITStepModel*>());
62 fReactionSet = nullptr; << 62
63 fpTrackContainer = nullptr; << 63 for(int i = 0; i < (int) G4ITType::size(); i++)
64 fpModelHandler = nullptr; << 64 {
65 fpActiveModelWithMinTimeStep = nullptr; << 65 fCurrentModel[i].assign(G4ITType::size(), 0);
66 fComputeTimeStep = false; << 66 }
67 fComputeReaction = false; << 67 fUserMinTimeStep = -1.;
>> 68 fTSTimeStep = DBL_MAX;
>> 69 fpTrackingManager = 0;
>> 70 fReactionSet = 0;
>> 71 fpTrackContainer = 0;
>> 72 fpModelHandler = 0;
>> 73 fComputeTimeStep = false;
>> 74 fComputeReaction = false;
>> 75 }
>> 76
>> 77 G4ITModelProcessor::~G4ITModelProcessor()
>> 78 {
>> 79 //dtor
>> 80 // if(fpModelHandler) delete fpModelHandler; deleted by G4Scheduler
>> 81 fCurrentModel.clear();
>> 82 fReactionInfo.clear();
>> 83 }
>> 84
>> 85 // Should not be used
>> 86 G4ITModelProcessor::G4ITModelProcessor(const G4ITModelProcessor& /*other*/)
>> 87 {
>> 88 //copy ctorr
>> 89 fpTrack = 0;
>> 90 fpModelHandler = 0;
>> 91 fpModel = 0;
>> 92 fInitialized = false;
>> 93 fpModelManager = 0;
>> 94 fUserMinTimeStep = -1.;
>> 95 fTSTimeStep = DBL_MAX;
>> 96 fpTrackingManager = 0;
>> 97 fReactionSet = 0;
>> 98 fpTrackContainer = 0;
>> 99 fComputeTimeStep = false;
>> 100 fComputeReaction = false;
68 } 101 }
69 102
70 G4ITModelProcessor::~G4ITModelProcessor() = de <<
71 <<
72 void G4ITModelProcessor::RegisterModel(double 103 void G4ITModelProcessor::RegisterModel(double time, G4VITStepModel* model)
73 { 104 {
74 fpModelHandler->RegisterModel(model, time) << 105 fpModelHandler->RegisterModel(model, time);
75 } 106 }
76 107
>> 108 // Should not be used
>> 109 G4ITModelProcessor& G4ITModelProcessor::operator=(const G4ITModelProcessor& rhs)
>> 110 {
>> 111 if(this == &rhs) return *this; // handle self assignment
>> 112 //assignment operator
>> 113 return *this;
>> 114 }
>> 115 //______________________________________________________________________________
77 void G4ITModelProcessor::Initialize() 116 void G4ITModelProcessor::Initialize()
78 { 117 {
79 fpModelHandler->Initialize(); << 118 fpModelHandler->Initialize();
80 fReactionSet = G4ITReactionSet::Instance() << 119 fReactionSet = G4ITReactionSet::Instance();
81 fpTrackContainer = G4ITTrackHolder::Instan << 120 fpTrackContainer = G4ITTrackHolder::Instance();
82 fInitialized = true; << 121 fInitialized = true;
83 fComputeTimeStep = false; << 122 fComputeTimeStep = false;
84 fComputeReaction = false; << 123 fComputeReaction = false;
85 if (fpModelHandler->GetTimeStepComputerFla << 124 if(fpModelHandler->GetTimeStepComputerFlag()) fComputeTimeStep = true;
86 { << 125 if(fpModelHandler->GetReactionProcessFlag()) fComputeReaction = true;
87 fComputeTimeStep = true; <<
88 } <<
89 if (fpModelHandler->GetReactionProcessFlag <<
90 { <<
91 fComputeReaction = true; <<
92 } <<
93 } 126 }
94 127
>> 128 //_________________________________________________________________________
>> 129
95 G4double G4ITModelProcessor::CalculateMinTimeS 130 G4double G4ITModelProcessor::CalculateMinTimeStep(G4double currentGlobalTime,
96 131 G4double definedMinTimeStep)
97 { 132 {
98 133
99 #if defined (DEBUG_MEM) && defined (DEBUG_MEM_ 134 #if defined (DEBUG_MEM) && defined (DEBUG_MEM_DETAILED_STEPPING)
100 MemStat mem_first, mem_second, mem_diff; << 135 MemStat mem_first, mem_second, mem_diff;
101 mem_first = MemoryUsage(); << 136 #endif
>> 137
>> 138 #if defined (DEBUG_MEM) && defined (DEBUG_MEM_DETAILED_STEPPING)
>> 139 mem_first = MemoryUsage();
102 #endif 140 #endif
103 141
104 fpActiveModelWithMinTimeStep = nullptr; << 142 fTSTimeStep = DBL_MAX;
105 fTSTimeStep = DBL_MAX; <<
106 143
107 InitializeStepper(currentGlobalTime, defin << 144 InitializeStepper(currentGlobalTime, definedMinTimeStep);
>> 145 // TODO
>> 146 // fpMasterModelProcessor -> InitializeStepper(fGlobalTime, fDefinedMinTimeStep) ;
108 147
109 #if defined (DEBUG_MEM) && defined (DEBUG_MEM_ 148 #if defined (DEBUG_MEM) && defined (DEBUG_MEM_DETAILED_STEPPING)
110 mem_second = MemoryUsage(); << 149 mem_second = MemoryUsage();
111 mem_diff = mem_second-mem_first; << 150 mem_diff = mem_second-mem_first;
112 G4cout << "\t || MEM || G4Scheduler::Calcu << 151 G4cout << "\t || MEM || G4Scheduler::CalculateMinTimeStep || After "
113 "computing fpModelProcessor -> InitializeS << 152 "computing fpModelProcessor -> InitializeStepper, diff is : "
114 << mem_diff << 153 << mem_diff
>> 154 << G4endl;
>> 155 #endif
>> 156
>> 157 // G4TrackList::iterator fpMainList_i = fpMainList->begin();
>> 158
>> 159 G4TrackManyList* mainList = fpTrackContainer->GetMainList();
>> 160 G4TrackManyList::iterator it = mainList->begin();
>> 161 G4TrackManyList::iterator end = mainList->end();
>> 162
>> 163 for (; it != end; ++it)
>> 164 {
>> 165 G4Track * track = *it;
>> 166
>> 167 if (track == 0)
>> 168 {
>> 169 G4ExceptionDescription exceptionDescription;
>> 170 exceptionDescription << "No track found.";
>> 171 G4Exception("G4Scheduler::CalculateMinStep", "ITScheduler006",
>> 172 FatalErrorInArgument, exceptionDescription);
>> 173 return 0; // makes coverity happy
>> 174 }
>> 175
>> 176 #ifdef DEBUG
>> 177 G4cout << "*_* " << GetIT(track)->GetName()
>> 178 << " ID: " << track->GetTrackID()
>> 179 << " at time : " << track->GetGlobalTime()
115 << G4endl; 180 << G4endl;
116 #endif 181 #endif
117 182
118 for (auto& pStepModel : fActiveModels) << 183 G4TrackStatus trackStatus = track->GetTrackStatus();
>> 184 if (trackStatus == fStopAndKill || trackStatus == fStopButAlive)
119 { 185 {
120 fTSTimeStep = << 186 continue;
121 pStepModel->GetTimeStepper()->Calc <<
122 currentGlobalTime, <<
123 definedMinTimeStep); <<
124 <<
125 fpActiveModelWithMinTimeStep = pStepMo <<
126 <<
127 if(fTSTimeStep == -1){ <<
128 fpActiveModelWithMinTimeStep->GetR <<
129 if(fReactionSet->Empty()) return D <<
130 const auto& fReactionSetInTime = f <<
131 fTSTimeStep = fReactionSetInTime.b <<
132 } <<
133 } 187 }
134 188
>> 189 // This Extract... was thought for MT mode at the track level
>> 190 //ExtractTimeStepperData(fpModelProcessor) ;
>> 191
>> 192 CalculateTimeStep(track, definedMinTimeStep);
>> 193
>> 194 // if MT mode at track level, this command should be displaced
>> 195 ExtractTimeStepperData();
>> 196 }
>> 197
135 #if defined (DEBUG_MEM) && defined (DEBUG_MEM_ 198 #if defined (DEBUG_MEM) && defined (DEBUG_MEM_DETAILED_STEPPING)
136 mem_second = MemoryUsage(); << 199 mem_second = MemoryUsage();
137 mem_diff = mem_second-mem_first; << 200 mem_diff = mem_second-mem_first;
138 G4cout << "\t || MEM || G4Scheduler::Calcu << 201 G4cout << "\t || MEM || G4Scheduler::CalculateMinTimeStep || "
139 "After looping on tracks, diff is : " << m << 202 "After looping on tracks, diff is : " << mem_diff << G4endl;
140 #endif 203 #endif
141 return fTSTimeStep; << 204 return fTSTimeStep;
>> 205 }
>> 206
>> 207 //_________________________________________________________________________
>> 208
>> 209 void G4ITModelProcessor::ExtractTimeStepperData()
>> 210 {
>> 211 if(fpTrack == 0)
>> 212 {
>> 213 CleanProcessor();
>> 214 return;
>> 215 }
>> 216
>> 217 const std::vector<std::vector<G4VITStepModel*> >* model = GetCurrentModel();
>> 218
>> 219 for(unsigned i = 0; i < model->size(); ++i)
>> 220 {
>> 221 for(unsigned j = 0; j < (*model)[i].size(); ++j)
>> 222 {
>> 223 G4VITStepModel* mod = (*model)[i][j];
>> 224
>> 225 if(mod == 0)
>> 226 {
>> 227 continue;
>> 228 }
>> 229
>> 230 G4VITTimeStepComputer* stepper(mod->GetTimeStepper());
>> 231
>> 232 G4double sampledMinTimeStep(stepper->GetSampledMinTimeStep());
>> 233 G4TrackVectorHandle reactants(stepper->GetReactants());
>> 234
>> 235 if(sampledMinTimeStep < fTSTimeStep)
>> 236 {
>> 237 /*
>> 238 // DEBUG SPECIAL CASE
>> 239 if(!reactants)
>> 240 {
>> 241 G4ExceptionDescription exceptionDescription ;
>> 242 CalculateMinTimeStep(); // => at least N (N = nb of tracks) loops exceptionDescription << "No reactants were found by the time stepper.";
>> 243 G4Exception("G4Scheduler::ExtractTimeStepperData","ITScheduler007",
>> 244 FatalErrorInArgument,exceptionDescription);
>> 245 continue ;
>> 246 }
>> 247 */
>> 248
>> 249 fTSTimeStep = sampledMinTimeStep;
>> 250 //fReactingTracks.clear();
>> 251
>> 252 fReactionSet->CleanAllReaction();
>> 253 if(bool(reactants))
>> 254 {
>> 255 // G4cout << "*** (1) G4Scheduler::ExtractTimeStepperData insert
>> 256 // reactants for " << GetIT(track)->GetName() << G4endl;
>> 257 // G4cout << "bool(reactants) = " << bool(reactants) << G4endl;
>> 258 // G4cout << reactants->size() << G4endl;
>> 259 // G4cout << GetIT(reactants->operator[](0))->GetName() << G4endl;
>> 260
>> 261 // fReactingTracks.insert(make_pair(track, reactants));
>> 262 fReactionSet->AddReactions(fTSTimeStep,
>> 263 const_cast<G4Track*>(fpTrack),
>> 264 reactants);
>> 265 stepper->ResetReactants();
>> 266 }
>> 267 }
>> 268 else if(fTSTimeStep == sampledMinTimeStep)
>> 269 {
>> 270 /*
>> 271 // DEBUG SPECIAL CASE
>> 272 if(!reactants)
>> 273 {
>> 274 G4ExceptionDescription exceptionDescription ;
>> 275 exceptionDescription << "No reactants were found by the time stepper.";
>> 276 G4Exception("G4Scheduler::ExtractTimeStepperData","ITScheduler008",
>> 277 FatalErrorInArgument,exceptionDescription);
>> 278 continue ;
>> 279 }
>> 280 */
>> 281 if(bool(reactants))
>> 282 {
>> 283 // G4cout << "*** (2) G4Scheduler::ExtractTimeStepperData insert
>> 284 // reactants for " << GetIT(track)->GetName() << G4endl;
>> 285 // G4cout << "bool(reactants) = " << bool(reactants) << G4endl;
>> 286 // G4cout << "trackA : " << GetIT(track)->GetName()
>> 287 // << " ("<< track->GetTrackID() << ")" << G4endl;
>> 288 // G4cout << reactants->size() << G4endl;
>> 289 // G4cout << GetIT(reactants->operator[](0))->GetName() << G4endl;
>> 290
>> 291 // fReactingTracks.insert(make_pair(track, reactants));
>> 292
>> 293 fReactionSet->AddReactions(fTSTimeStep,
>> 294 const_cast<G4Track*>(fpTrack),
>> 295 reactants);
>> 296 stepper->ResetReactants();
>> 297 }
>> 298 }
>> 299 else
>> 300 {
>> 301 if(bool(reactants))
>> 302 {
>> 303 stepper->ResetReactants();
>> 304 }
>> 305 }
>> 306 }
>> 307 }
>> 308
>> 309 CleanProcessor();
142 } 310 }
143 311
144 //____________________________________________ 312 //______________________________________________________________________________
145 313
146 void G4ITModelProcessor::InitializeStepper(G4d 314 void G4ITModelProcessor::InitializeStepper(G4double currentGlobalTime,
147 G4d 315 G4double userMinTime)
148 { 316 {
149 G4VITTimeStepComputer::SetTimes(currentGlo << 317 // G4cout << "G4ITModelProcessor::InitializeStepper" << G4endl;
>> 318 if(fpModelHandler == 0)
>> 319 {
>> 320 G4ExceptionDescription exceptionDescription;
>> 321 exceptionDescription
>> 322 << "No G4ITModelHandler was passed to the modelProcessor.";
>> 323 G4Exception("G4ITModelProcessor::InitializeStepper",
>> 324 "ITModelProcessor002",
>> 325 FatalErrorInArgument,
>> 326 exceptionDescription);
>> 327 }
>> 328 const std::vector<std::vector<G4ITModelManager*> >* modelManager =
>> 329 fpModelHandler->GetAllModelManager();
>> 330
>> 331 if(modelManager == 0)
>> 332 {
>> 333 G4ExceptionDescription exceptionDescription;
>> 334 exceptionDescription
>> 335 << "No G4ITModelManager was register to G4ITModelHandler.";
>> 336 G4Exception("G4ITModelProcessor::InitializeStepper",
>> 337 "ITModelProcessor003",
>> 338 FatalErrorInArgument,
>> 339 exceptionDescription);
>> 340 }
>> 341
>> 342 int nbModels1 = modelManager->size();
>> 343
>> 344 G4VITTimeStepComputer::SetTimes(currentGlobalTime, userMinTime);
>> 345
>> 346 // TODO !!!
>> 347 // if( nbModels1 != 1 || (nbModels1 == 1 && !fpModelManager) )
>> 348 {
>> 349 int nbModels2 = -1;
>> 350 G4VITStepModel* model = 0;
>> 351 G4ITModelManager* modman = 0;
>> 352
>> 353 for(int i = 0; i < nbModels1; i++)
>> 354 {
>> 355 nbModels2 = (*modelManager)[i].size();
>> 356
>> 357 for(int j = 0; j <= i; j++)
>> 358 {
>> 359 modman = (*modelManager)[i][j];
>> 360
>> 361 if(modman == 0) continue;
>> 362
>> 363 model = modman->GetModel(currentGlobalTime);
>> 364 G4VITTimeStepComputer* stepper = model->GetTimeStepper();
150 365
151 #if defined (DEBUG_MEM) 366 #if defined (DEBUG_MEM)
152 MemStat mem_first, mem_second, mem_diff; << 367 MemStat mem_first, mem_second, mem_diff;
153 mem_first = MemoryUsage(); <<
154 #endif 368 #endif
155 369
156 fActiveModels = fpModelHandler->GetActiveM << 370 #if defined (DEBUG_MEM)
>> 371 mem_first = MemoryUsage();
>> 372 #endif
157 373
158 for (auto& pModel : fActiveModels) << 374 // stepper->PrepareForAllProcessors() ;
159 { << 375 stepper->Prepare();
160 pModel->PrepareNewTimeStep(); <<
161 } <<
162 376
163 #if defined (DEBUG_MEM) 377 #if defined (DEBUG_MEM)
164 mem_second = MemoryUsage(); << 378 mem_second = MemoryUsage();
165 mem_diff = mem_second-mem_first; << 379 mem_diff = mem_second-mem_first;
166 G4cout << "\t || MEM || G4ITModelP << 380 G4cout << "\t || MEM || G4ITModelProcessor::InitializeStepper || After computing stepper -> Prepare(), diff is : " << mem_diff << G4endl;
167 #endif 381 #endif
>> 382 fCurrentModel[i][j] = model;
>> 383 }
>> 384 }
>> 385
>> 386 if(nbModels1 == 1 && nbModels2 == 1)
>> 387 {
>> 388 fpModelManager = modman;
>> 389 fpModel = model;
>> 390 }
>> 391 else fpModel = 0;
>> 392 }
>> 393 }
>> 394
>> 395 //______________________________________________________________________________
>> 396 void G4ITModelProcessor::CalculateTimeStep(const G4Track* track,
>> 397 const G4double userMinTimeStep)
>> 398 {
>> 399 // G4cout << "G4ITModelProcessor::CalculateStep" << G4endl;
>> 400 CleanProcessor();
>> 401 if(track == 0)
>> 402 {
>> 403 G4ExceptionDescription exceptionDescription;
>> 404 exceptionDescription << "No track was passed to the method (track == 0).";
>> 405 G4Exception("G4ITModelProcessor::CalculateStep",
>> 406 "ITModelProcessor004",
>> 407 FatalErrorInArgument,
>> 408 exceptionDescription);
>> 409 }
>> 410 SetTrack(track);
>> 411 fUserMinTimeStep = userMinTimeStep;
>> 412
>> 413 DoCalculateStep();
>> 414 }
168 415
>> 416 //______________________________________________________________________________
>> 417 void G4ITModelProcessor::DoCalculateStep()
>> 418 {
>> 419 if(fpModel) // ie only one model has been declared and will be used
>> 420 {
>> 421 fpModel->GetTimeStepper()->CalculateStep(*fpTrack, fUserMinTimeStep);
>> 422 }
>> 423 else // ie many models have been declared and will be used
>> 424 {
>> 425 std::vector<G4VITStepModel*>& model = fCurrentModel[GetIT(fpTrack)
>> 426 ->GetITType()];
>> 427
>> 428 for(int i = 0; i < (int) model.size(); i++)
>> 429 {
>> 430 if(model[i] == 0) continue;
>> 431 model[i]->GetTimeStepper()->CalculateStep(*fpTrack, fUserMinTimeStep);
>> 432 }
>> 433 }
169 } 434 }
170 435
171 //____________________________________________ 436 //_________________________________________________________________________
172 437
173 void G4ITModelProcessor::ComputeTrackReaction( 438 void G4ITModelProcessor::ComputeTrackReaction(G4ITStepStatus fITStepStatus,
174 439 G4double fGlobalTime,
175 440 G4double currentTimeStep,
176 << 441 G4double previousTimeStep,
177 442 G4bool reachedUserTimeLimit,
178 443 G4double fTimeTolerance,
179 444 G4UserTimeStepAction* fpUserTimeStepAction,
180 445 G4int
181 #ifdef G4VERBOSE 446 #ifdef G4VERBOSE
182 fVerbose << 447 fVerbose
183 #endif 448 #endif
184 ) << 449 )
185 { 450 {
186 if (fReactionSet->Empty()) << 451 // if (fReactingTracks.empty())
187 { << 452 if(fReactionSet->Empty())
188 return; << 453 {
189 } << 454 return;
>> 455 }
>> 456
>> 457 if(fITStepStatus == eCollisionBetweenTracks)
>> 458 // if(fInteractionStep == false)
>> 459 {
>> 460 // TODO
>> 461 FindReaction(fReactionSet,
>> 462 currentTimeStep,
>> 463 previousTimeStep,
>> 464 reachedUserTimeLimit);
>> 465 // TODO
>> 466 // A ne faire uniquement si le temps choisis est celui calculé par le time stepper
>> 467 // Sinon utiliser quelque chose comme : fModelProcessor->FindReaction(&fMainList);
>> 468
>> 469 std::vector<G4ITReactionChange*> * reactionInfo_v = GetReactionInfo();
>> 470 std::vector<G4ITReactionChange*>::iterator reactionInfo_i = reactionInfo_v
>> 471 ->begin();
190 472
191 if (fITStepStatus == eCollisionBetweenTrac << 473 for(; reactionInfo_i != reactionInfo_v->end(); ++reactionInfo_i)
192 { 474 {
193 G4VITReactionProcess* pReactionProcess << 475 G4ITReactionChange* changes = (*reactionInfo_i);
194 fReactionInfo = pReactionProcess->Find << 476 G4Track* trackA = const_cast<G4Track*>(changes->GetTrackA());
195 currentTimeStep, << 477 G4Track* trackB = const_cast<G4Track*>(changes->GetTrackB());
196 fGlobalTime, << 478
197 reachedUserTimeLimit); << 479 if(trackA == 0 || trackB == 0 || trackA->GetTrackStatus() == fStopAndKill
198 << 480 || trackB->GetTrackStatus() == fStopAndKill) continue;
199 // TODO << 481
200 // A ne faire uniquement si le temps c << 482 G4int nbSecondaries = changes->GetNumberOfSecondaries();
201 // Sinon utiliser quelque chose comme << 483 const vector<G4Track*>* productsVector = changes->GetfSecondary();
202 << 484
203 for (auto& pChanges : fReactionInfo) << 485 if(fpUserTimeStepAction)
204 { << 486 {
205 auto pTrackA = const_cast<G4Track* << 487 fpUserTimeStepAction->UserReactionAction(*trackA,
206 auto pTrackB = const_cast<G4Track* << 488 *trackB,
207 << 489 productsVector);
208 if (pTrackA == nullptr << 490 }
209 || pTrackB == nullptr <<
210 || pTrackA->GetTrackStatus() = <<
211 || pTrackB->GetTrackStatus() = <<
212 { <<
213 continue; <<
214 } <<
215 <<
216 G4int nbSecondaries = pChanges->Ge <<
217 const std::vector<G4Track*>* produ <<
218 <<
219 if (fpUserTimeStepAction != nullpt <<
220 { <<
221 fpUserTimeStepAction->UserReac <<
222 <<
223 <<
224 } <<
225 491
226 #ifdef G4VERBOSE 492 #ifdef G4VERBOSE
227 if (fVerbose != 0) << 493 if(fVerbose)
228 { << 494 {
229 G4cout << "At time : " << std: << 495 G4cout << "At time : " << setw(7) << G4BestUnit(fGlobalTime, "Time")
230 << " Reaction : " << Ge << 496 << " Reaction : " << GetIT(trackA)->GetName() << " ("
231 << pTrackA->GetTrackID( << 497 << trackA->GetTrackID() << ") + " << GetIT(trackB)->GetName() << " ("
232 << pTrackB->GetTrackID( << 498 << trackB->GetTrackID() << ") -> ";
233 } << 499 }
234 #endif 500 #endif
235 501
236 if (nbSecondaries > 0) << 502 if(nbSecondaries > 0)
237 { << 503 {
238 for (int i = 0; i < nbSecondar << 504 for(int i = 0; i < nbSecondaries; ++i)
239 { << 505 {
240 #ifdef G4VERBOSE 506 #ifdef G4VERBOSE
241 if ((fVerbose != 0) && i ! << 507 if(fVerbose && i != 0) G4cout << " + ";
242 { << 508 #endif
243 G4cout << " + "; << 509
244 } << 510 G4Track* secondary = (*productsVector)[i]; //changes->GetSecondary(i);
245 #endif << 511 fpTrackContainer->_PushTrack(secondary);
246 << 512 GetIT(secondary)->SetParentID(trackA->GetTrackID(),
247 G4Track* secondary = (*pro << 513 trackB->GetTrackID());
248 // fpTrackContainer->_PushT << 514
249 GetIT(secondary)->SetParen << 515 if(secondary->GetGlobalTime() - fGlobalTime > fTimeTolerance)
250 << 516 {
251 << 517 G4ExceptionDescription exceptionDescription;
252 if (secondary->GetGlobalTi << 518 exceptionDescription << "The time of the secondary should not be bigger than the"
253 { << 519 " current global time."
254 G4ExceptionDescription << 520 << " This may cause synchronization problem. If the process you"
255 exceptionDescription < << 521 " are using required "
256 << 522 << "such feature please contact the developpers." << G4endl<< "The global time in the step manager : "
257 < << 523 << G4BestUnit(fGlobalTime,"Time")
258 << 524 << G4endl
259 < << 525 << "The global time of the track : "
260 < << 526 << G4BestUnit(secondary->GetGlobalTime(),"Time")
261 < << 527 << G4endl;
262 < << 528
263 < << 529 G4Exception("G4Scheduler::ComputeInteractionBetweenTracks",
264 < << 530 "ITScheduler010",
265 < << 531 FatalErrorInArgument,
266 << 532 exceptionDescription);
267 G4Exception("G4Schedul << 533 }
268 "ITSchedul <<
269 FatalError <<
270 exceptionD <<
271 } <<
272 534
273 #ifdef G4VERBOSE 535 #ifdef G4VERBOSE
274 if (fVerbose != 0) << 536 if(fVerbose)
275 { << 537 G4cout << GetIT(secondary)->GetName() << " ("
276 G4cout << GetIT(second << 538 << secondary->GetTrackID() << ")";
277 << secondary->G << 539 #endif
278 } << 540 }
279 #endif << 541 }
280 } << 542 else
281 } << 543 {
282 else <<
283 { <<
284 #ifdef G4VERBOSE 544 #ifdef G4VERBOSE
285 if (fVerbose != 0) << 545 if(fVerbose) G4cout << "No product";
286 { <<
287 G4cout << "No product"; <<
288 } <<
289 #endif 546 #endif
290 } << 547 }
291 #ifdef G4VERBOSE 548 #ifdef G4VERBOSE
292 if (fVerbose != 0) << 549 if(fVerbose) G4cout << G4endl;
293 { << 550 #endif
294 G4cout << G4endl; << 551 if(trackA->GetTrackID() == 0 || trackB->GetTrackID() == 0)
295 } << 552 {
296 #endif << 553 G4Track* track = 0;
297 if (pTrackA->GetTrackID() == 0 || << 554 if(trackA->GetTrackID() == 0) track = trackA;
298 { << 555 else track = trackB;
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 } <<
335 <<
336 if (pChanges->WereParentsKilled()) <<
337 { <<
338 pTrackA->SetTrackStatus(fStopA <<
339 pTrackB->SetTrackStatus(fStopA <<
340 <<
341 fpTrackingManager->EndTracking <<
342 fpTrackingManager->EndTracking <<
343 } <<
344 556
345 pChanges.reset(nullptr); << 557 G4ExceptionDescription exceptionDescription;
>> 558 exceptionDescription
>> 559 << "The problem was found for the reaction between tracks :"
>> 560 << trackA->GetParticleDefinition()->GetParticleName() << " ("
>> 561 << trackA->GetTrackID() << ") & "
>> 562 << trackB->GetParticleDefinition()->GetParticleName() << " ("
>> 563 << trackB->GetTrackID() << "). \n";
>> 564
>> 565 if(track->GetStep() == 0)
>> 566 {
>> 567 exceptionDescription << "Also no step was found"
>> 568 << " ie track->GetStep() == 0 \n";
346 } 569 }
347 570
348 fReactionInfo.clear(); << 571 exceptionDescription << "Parent ID of trackA : "
>> 572 << trackA->GetParentID() << "\n";
>> 573 exceptionDescription << "Parent ID of trackB : "
>> 574 << trackB->GetParentID() << "\n";
>> 575
>> 576 exceptionDescription
>> 577 << "The ID of one of the reaction track was not setup.";
>> 578 G4Exception("G4Scheduler::ComputeInteractionBetweenTracks",
>> 579 "ITScheduler011",
>> 580 FatalErrorInArgument,
>> 581 exceptionDescription);
>> 582 }
>> 583
>> 584 if(changes->WereParentsKilled())
>> 585 {
>> 586 // DEBUG
>> 587 // G4cout << "Erasing tracks : "
>> 588 // << "trackA at time : " << G4BestUnit(trackA->GetGlobalTime() , "Time")
>> 589 // << "\t trackB at time : "<< G4BestUnit(trackB->GetGlobalTime(), "Time")
>> 590 // << "\t GlobalTime : " << G4BestUnit(fGlobalTime, "Time")
>> 591 // << G4endl;
>> 592
>> 593 trackA->SetTrackStatus(fStopAndKill);
>> 594 trackB->SetTrackStatus(fStopAndKill);
>> 595
>> 596 // G4TrackList::Pop(trackA);
>> 597 // G4TrackList::Pop(trackB);
>> 598 fpTrackingManager->EndTracking(trackA);
>> 599 fpTrackingManager->EndTracking(trackB);
>> 600 }
>> 601
>> 602 delete changes;
349 } 603 }
350 604
351 // fReactionSet->CleanAllReaction(); << 605 reactionInfo_v->clear(); // never null because point to a non-pointer member of ModelProcessor
>> 606 }
>> 607 // DEBUG
>> 608 //else
>> 609 //{
>> 610 // G4cout << "fInteractionStep == true" << G4endl ;
>> 611 //}
352 612
353 fpTrackContainer->MergeSecondariesWithMain << 613 // fReactingTracks.clear();
354 fpTrackContainer->KillTracks(); << 614 fReactionSet->CleanAllReaction();
355 } <<
356 615
357 void G4ITModelProcessor::SetTrack(const G4Trac << 616 fpTrackContainer->MergeSecondariesWithMainList();
358 { << 617 fpTrackContainer->KillTracks();
359 fpTrack = track; <<
360 } 618 }
361 619
362 void G4ITModelProcessor::SetModelHandler(G4ITM << 620 //______________________________________________________________________________
363 { << 621 void G4ITModelProcessor::FindReaction(G4ITReactionSet* reactionSet,
364 if (fInitialized) << 622 const double currentStepTime,
>> 623 const double previousStepTime,
>> 624 const bool reachedUserStepTimeLimit)
>> 625 {
>> 626 // DEBUG
>> 627 // G4cout << "G4ITReactionManager::FindReaction" << G4endl;
>> 628 //if (tracks == 0) return;
>> 629 if(reactionSet == 0) return;
>> 630 if(fpModelHandler->GetAllModelManager()->empty()) return;
>> 631
>> 632 G4ITReactionPerTrackMap& reactionPerTrackMap = reactionSet->GetReactionMap();
>> 633
>> 634 std::map<G4Track*, G4ITReactionPerTrackPtr, compTrackPerID>::iterator tracks_i =
>> 635 reactionPerTrackMap.begin();
>> 636
>> 637 //std::map<G4Track*, G4TrackVectorHandle, compTrackPerID>::iterator tracks_i = tracks->begin();
>> 638
>> 639 // G4cout << "G4ITModelProcessor::FindReaction at step :" << G4ITTimeStepper::Instance()->GetNbSteps() << G4endl;
>> 640
>> 641 // for (tracks_i = tracks->begin(); tracks_i != tracks->end(); tracks_i++)
>> 642 for(tracks_i = reactionPerTrackMap.begin();
>> 643 tracks_i != reactionPerTrackMap.end();
>> 644 tracks_i = reactionPerTrackMap.begin())
>> 645 {
>> 646 //G4cout << "here" << G4endl;
>> 647 G4Track* trackA = tracks_i->first;
>> 648 if(trackA->GetTrackStatus() == fStopAndKill)
365 { 649 {
366 G4ExceptionDescription exceptionDescri << 650 //G4cout << "continue 1" << G4endl;
367 exceptionDescription << 651 continue;
368 << "You are trying to set a new mo <<
369 G4Exception("G4ITModelProcessor::SetMo <<
370 FatalErrorInArgument, exce <<
371 } 652 }
372 fpModelHandler = pModelHandler; <<
373 } <<
374 653
375 void G4ITModelProcessor::CleanProcessor() << 654 G4ITReactionPerTrackPtr reactionPerTrack = tracks_i->second;
376 { << 655 G4ITReactionList& reactionList = reactionPerTrack->GetReactionList();
377 fpTrack = nullptr; <<
378 } <<
379 656
380 bool G4ITModelProcessor::GetComputeTimeStep() << 657 G4IT* ITA = GetIT(trackA);
381 { << 658 G4ITType ITypeA = ITA->GetITType();
382 return fComputeTimeStep; <<
383 } <<
384 659
385 const G4Track* G4ITModelProcessor::GetTrack() << 660 const std::vector<G4VITStepModel*> model = fCurrentModel[ITypeA];
386 { <<
387 return fpTrack; <<
388 } <<
389 661
390 void G4ITModelProcessor::SetTrackingManager(G4 << 662 G4Track* trackB = 0;
391 { << 663 G4ITType ITypeB(-1);
392 fpTrackingManager = pTrackingManager; << 664 G4VITReactionProcess* process = 0;
393 } << 665 G4ITReactionChange* changes = 0;
>> 666
>> 667 assert(reactionList.begin() != reactionList.end());
394 668
>> 669 for(G4ITReactionList::iterator it = reactionList.begin();
>> 670 it != reactionList.end(); it = reactionList.begin())
>> 671 {
>> 672 G4ITReactionPtr reaction(*it);
>> 673 trackB = reaction->GetReactant(trackA);
>> 674 if(trackB->GetTrackStatus() == fStopAndKill)
>> 675 {
>> 676 //G4cout << "continue 2" << G4endl;
>> 677 continue;
>> 678 }
>> 679
>> 680 if(trackB == trackA)
>> 681 {
>> 682 G4ExceptionDescription exceptionDescription;
>> 683 exceptionDescription
>> 684 << "The IT reaction process sent back a reaction between trackA and trackB. ";
>> 685 exceptionDescription << "The problem is trackA == trackB";
>> 686 G4Exception("G4ITModelProcessor::FindReaction",
>> 687 "ITModelProcessor005",
>> 688 FatalErrorInArgument,
>> 689 exceptionDescription);
>> 690 }
>> 691
>> 692 G4IT* ITB = GetIT(trackB);
>> 693 G4ITType ITypeBtmp = ITB->GetITType();
>> 694
>> 695 if(ITypeB != ITypeBtmp)
>> 696 {
>> 697 ITypeB = ITypeBtmp;
>> 698
>> 699 if(model[ITypeB]) process = model[ITypeB]->GetReactionProcess();
>> 700 }
>> 701
>> 702 reactionSet->SelectThisReaction(reaction);
>> 703
>> 704 if(process && process->TestReactibility(*trackA,
>> 705 *trackB,
>> 706 currentStepTime,
>> 707 previousStepTime,
>> 708 reachedUserStepTimeLimit))
>> 709 {
>> 710 changes = process->MakeReaction(*trackA, *trackB);
>> 711 }
>> 712
>> 713 if(changes)
>> 714 {
>> 715 fReactionInfo.push_back(changes);
>> 716
>> 717 // G4cout << "pushing reaction for trackA (" << trackA->GetTrackID() << ") and trackB ("
>> 718 // << trackB->GetTrackID() << ")" << G4endl;
>> 719 //
>> 720 // G4cout << "nb of secondaries : " << changes->GetNumberOfSecondaries() << G4endl;
>> 721 // G4cout << "with track 0 = " << changes->GetSecondary(0) << G4endl;
>> 722
>> 723 process->ResetChanges();
>> 724 changes = 0;
>> 725
>> 726 break;
>> 727 }
>> 728 }
>> 729 }
>> 730
>> 731 //assert(G4ITReaction::gAll->empty() == true);
>> 732 }
395 733