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