Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/electromagnetic/dna/management/src/G4ITModelProcessor.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

  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