Geant4 Cross Reference |
1 // 2 // ******************************************************************** 3 // * License and Disclaimer * 4 // * * 5 // * The Geant4 software is copyright of the Copyright Holders of * 6 // * the Geant4 Collaboration. It is provided under the terms and * 7 // * conditions of the Geant4 Software License, included in the file * 8 // * LICENSE and available at http://cern.ch/geant4/license . These * 9 // * include a list of copyright holders. * 10 // * * 11 // * Neither the authors of this software system, nor their employing * 12 // * institutes,nor the agencies providing financial support for this * 13 // * work make any representation or warranty, express or implied, * 14 // * regarding this software system or assume any liability for its * 15 // * use. Please see the license in the file LICENSE and URL above * 16 // * for the full disclaimer and the limitation of liability. * 17 // * * 18 // * This code implementation is the result of the scientific and * 19 // * technical work of the GEANT4 collaboration. * 20 // * By using, copying, modifying or distributing the software (or * 21 // * any work based on the software) you agree to acknowledge its * 22 // * use in resulting scientific publications, and indicate your * 23 // * acceptance of all terms of the Geant4 Software license. * 24 // ******************************************************************** 25 // 26 // 27 // Author: Mathieu Karamitros (kara (AT) cenbg . in2p3 . fr) 28 // 29 // History: 30 // ----------- 31 // 10 Oct 2011 M.Karamitros created 32 // 33 // ------------------------------------------------------------------- 34 35 #include "G4ITModelProcessor.hh" 36 #include "G4VITTimeStepComputer.hh" 37 #include "G4VITReactionProcess.hh" 38 #include "G4ITReaction.hh" 39 #include "G4ITTrackHolder.hh" 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 55 G4ITModelProcessor::G4ITModelProcessor() 56 { 57 fpTrack = nullptr; 58 fInitialized = false; 59 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 } 69 70 G4ITModelProcessor::~G4ITModelProcessor() = default; 71 72 void G4ITModelProcessor::RegisterModel(double time, G4VITStepModel* model) 73 { 74 fpModelHandler->RegisterModel(model, time); 75 } 76 77 void G4ITModelProcessor::Initialize() 78 { 79 fpModelHandler->Initialize(); 80 fReactionSet = G4ITReactionSet::Instance(); 81 fpTrackContainer = G4ITTrackHolder::Instance(); 82 fInitialized = true; 83 fComputeTimeStep = false; 84 fComputeReaction = false; 85 if (fpModelHandler->GetTimeStepComputerFlag()) 86 { 87 fComputeTimeStep = true; 88 } 89 if (fpModelHandler->GetReactionProcessFlag()) 90 { 91 fComputeReaction = true; 92 } 93 } 94 95 G4double G4ITModelProcessor::CalculateMinTimeStep(G4double currentGlobalTime, 96 G4double definedMinTimeStep) 97 { 98 99 #if defined (DEBUG_MEM) && defined (DEBUG_MEM_DETAILED_STEPPING) 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, definedMinTimeStep); 108 109 #if defined (DEBUG_MEM) && defined (DEBUG_MEM_DETAILED_STEPPING) 110 mem_second = MemoryUsage(); 111 mem_diff = mem_second-mem_first; 112 G4cout << "\t || MEM || G4Scheduler::CalculateMinTimeStep || After " 113 "computing fpModelProcessor -> InitializeStepper, diff is : " 114 << mem_diff 115 << G4endl; 116 #endif 117 118 for (auto& pStepModel : fActiveModels) 119 { 120 fTSTimeStep = 121 pStepModel->GetTimeStepper()->CalculateMinTimeStep( 122 currentGlobalTime, 123 definedMinTimeStep); 124 125 fpActiveModelWithMinTimeStep = pStepModel; 126 127 if(fTSTimeStep == -1){ 128 fpActiveModelWithMinTimeStep->GetReactionProcess()->Initialize(); 129 if(fReactionSet->Empty()) return DBL_MAX; 130 const auto& fReactionSetInTime = fReactionSet->GetReactionsPerTime(); 131 fTSTimeStep = fReactionSetInTime.begin()->get()->GetTime() - currentGlobalTime; 132 } 133 } 134 135 #if defined (DEBUG_MEM) && defined (DEBUG_MEM_DETAILED_STEPPING) 136 mem_second = MemoryUsage(); 137 mem_diff = mem_second-mem_first; 138 G4cout << "\t || MEM || G4Scheduler::CalculateMinTimeStep || " 139 "After looping on tracks, diff is : " << mem_diff << G4endl; 140 #endif 141 return fTSTimeStep; 142 } 143 144 //______________________________________________________________________________ 145 146 void G4ITModelProcessor::InitializeStepper(G4double currentGlobalTime, 147 G4double userMinTime) 148 { 149 G4VITTimeStepComputer::SetTimes(currentGlobalTime, userMinTime); 150 151 #if defined (DEBUG_MEM) 152 MemStat mem_first, mem_second, mem_diff; 153 mem_first = MemoryUsage(); 154 #endif 155 156 fActiveModels = fpModelHandler->GetActiveModels(currentGlobalTime); 157 158 for (auto& pModel : fActiveModels) 159 { 160 pModel->PrepareNewTimeStep(); 161 } 162 163 #if defined (DEBUG_MEM) 164 mem_second = MemoryUsage(); 165 mem_diff = mem_second-mem_first; 166 G4cout << "\t || MEM || G4ITModelProcessor::InitializeStepper || After computing stepper -> Prepare(), diff is : " << mem_diff << G4endl; 167 #endif 168 169 } 170 171 //_________________________________________________________________________ 172 173 void G4ITModelProcessor::ComputeTrackReaction(G4ITStepStatus fITStepStatus, 174 G4double fGlobalTime, 175 G4double currentTimeStep, 176 G4double /*previousTimeStep*/, 177 G4bool reachedUserTimeLimit, 178 G4double fTimeTolerance, 179 G4UserTimeStepAction* fpUserTimeStepAction, 180 G4int 181 #ifdef G4VERBOSE 182 fVerbose 183 #endif 184 ) 185 { 186 if (fReactionSet->Empty()) 187 { 188 return; 189 } 190 191 if (fITStepStatus == eCollisionBetweenTracks) 192 { 193 G4VITReactionProcess* pReactionProcess = fpActiveModelWithMinTimeStep->GetReactionProcess(); 194 fReactionInfo = pReactionProcess->FindReaction(fReactionSet, 195 currentTimeStep, 196 fGlobalTime, 197 reachedUserTimeLimit); 198 199 // TODO 200 // A ne faire uniquement si le temps choisis est celui calculé par le time stepper 201 // Sinon utiliser quelque chose comme : fModelProcessor->FindReaction(&fMainList); 202 203 for (auto& pChanges : fReactionInfo) 204 { 205 auto pTrackA = const_cast<G4Track*>(pChanges->GetTrackA()); 206 auto pTrackB = const_cast<G4Track*>(pChanges->GetTrackB()); 207 208 if (pTrackA == nullptr 209 || pTrackB == nullptr 210 || pTrackA->GetTrackStatus() == fStopAndKill 211 || pTrackB->GetTrackStatus() == fStopAndKill) 212 { 213 continue; 214 } 215 216 G4int nbSecondaries = pChanges->GetNumberOfSecondaries(); 217 const std::vector<G4Track*>* productsVector = pChanges->GetfSecondary(); 218 219 if (fpUserTimeStepAction != nullptr) 220 { 221 fpUserTimeStepAction->UserReactionAction(*pTrackA, 222 *pTrackB, 223 productsVector); 224 } 225 226 #ifdef G4VERBOSE 227 if (fVerbose != 0) 228 { 229 G4cout << "At time : " << std::setw(7) << G4BestUnit(fGlobalTime, "Time") 230 << " Reaction : " << GetIT(pTrackA)->GetName() << " (" 231 << pTrackA->GetTrackID() << ") + " << GetIT(pTrackB)->GetName() << " (" 232 << pTrackB->GetTrackID() << ") -> "; 233 } 234 #endif 235 236 if (nbSecondaries > 0) 237 { 238 for (int i = 0; i < nbSecondaries; ++i) 239 { 240 #ifdef G4VERBOSE 241 if ((fVerbose != 0) && i != 0) 242 { 243 G4cout << " + "; 244 } 245 #endif 246 247 G4Track* secondary = (*productsVector)[i]; //changes->GetSecondary(i); 248 // fpTrackContainer->_PushTrack(secondary); 249 GetIT(secondary)->SetParentID(pTrackA->GetTrackID(), 250 pTrackB->GetTrackID()); 251 252 if (secondary->GetGlobalTime() - fGlobalTime > fTimeTolerance) 253 { 254 G4ExceptionDescription exceptionDescription; 255 exceptionDescription << "The time of the secondary should not be bigger than the" 256 " current global time." 257 << " This may cause synchronization problem. If the process you" 258 " are using required " 259 << "such feature please contact the developers." << G4endl 260 << "The global time in the step manager : " 261 << G4BestUnit(fGlobalTime, "Time") 262 << G4endl 263 << "The global time of the track : " 264 << G4BestUnit(secondary->GetGlobalTime(), "Time") 265 << G4endl; 266 267 G4Exception("G4Scheduler::ComputeInteractionBetweenTracks", 268 "ITScheduler010", 269 FatalErrorInArgument, 270 exceptionDescription); 271 } 272 273 #ifdef G4VERBOSE 274 if (fVerbose != 0) 275 { 276 G4cout << GetIT(secondary)->GetName() << " (" 277 << secondary->GetTrackID() << ")"; 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 || pTrackB->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 exceptionDescription; 310 exceptionDescription 311 << "The problem was found for the reaction between tracks :" 312 << pTrackA->GetParticleDefinition()->GetParticleName() << " (" 313 << pTrackA->GetTrackID() << ") & " 314 << pTrackB->GetParticleDefinition()->GetParticleName() << " (" 315 << pTrackB->GetTrackID() << "). \n"; 316 317 if (pTrack->GetStep() == nullptr) 318 { 319 exceptionDescription << "Also no step was found" 320 << " ie track->GetStep() == 0 \n"; 321 } 322 323 exceptionDescription << "Parent ID of trackA : " 324 << pTrackA->GetParentID() << "\n"; 325 exceptionDescription << "Parent ID of trackB : " 326 << pTrackB->GetParentID() << "\n"; 327 328 exceptionDescription 329 << "The ID of one of the reaction track was not setup."; 330 G4Exception("G4Scheduler::ComputeInteractionBetweenTracks", 331 "ITScheduler011", 332 FatalErrorInArgument, 333 exceptionDescription); 334 } 335 336 if (pChanges->WereParentsKilled()) 337 { 338 pTrackA->SetTrackStatus(fStopAndKill); 339 pTrackB->SetTrackStatus(fStopAndKill); 340 341 fpTrackingManager->EndTracking(pTrackA); 342 fpTrackingManager->EndTracking(pTrackB); 343 } 344 345 pChanges.reset(nullptr); 346 } 347 348 fReactionInfo.clear(); 349 } 350 351 // fReactionSet->CleanAllReaction(); 352 353 fpTrackContainer->MergeSecondariesWithMainList(); 354 fpTrackContainer->KillTracks(); 355 } 356 357 void G4ITModelProcessor::SetTrack(const G4Track* track) 358 { 359 fpTrack = track; 360 } 361 362 void G4ITModelProcessor::SetModelHandler(G4ITModelHandler* pModelHandler) 363 { 364 if (fInitialized) 365 { 366 G4ExceptionDescription exceptionDescription; 367 exceptionDescription 368 << "You are trying to set a new model while the model processor has alreaday be initialized"; 369 G4Exception("G4ITModelProcessor::SetModelHandler", "ITModelProcessor001", 370 FatalErrorInArgument, exceptionDescription); 371 } 372 fpModelHandler = pModelHandler; 373 } 374 375 void G4ITModelProcessor::CleanProcessor() 376 { 377 fpTrack = nullptr; 378 } 379 380 bool G4ITModelProcessor::GetComputeTimeStep() const 381 { 382 return fComputeTimeStep; 383 } 384 385 const G4Track* G4ITModelProcessor::GetTrack() const 386 { 387 return fpTrack; 388 } 389 390 void G4ITModelProcessor::SetTrackingManager(G4ITTrackingManager* pTrackingManager) 391 { 392 fpTrackingManager = pTrackingManager; 393 } 394 395