| 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 "G4ITStepProcessor.hh"
36 #include "G4LossTableManager.hh"
37 #include "G4EnergyLossTables.hh"
38 #include "G4ProductionCuts.hh"
39 #include "G4ProductionCutsTable.hh"
40 #include "G4VITProcess.hh"
41 #include "G4TrackingInformation.hh"
42 #include "G4IT.hh"
43 #include "G4ITTrackingManager.hh"
44 #include "G4ITTransportation.hh"
45
46 #include "G4ITNavigator.hh" // Include from 'geometry'
47
48 #include "G4ITSteppingVerbose.hh"
49 #include "G4VITSteppingVerbose.hh"
50
51 #include "G4ITTrackHolder.hh"
52 #include "G4ITReaction.hh"
53
54
55 //#define DEBUG_MEM 1
56
57 #ifdef DEBUG_MEM
58 #include "G4MemStat.hh"
59 using namespace G4MemStat;
60 using G4MemStat::MemStat;
61 #endif
62
63 void G4ITStepProcessor::DealWithSecondaries(G4int& counter)
64 {
65 // Now Store the secondaries from ParticleChange to SecondaryList
66 G4Track* tempSecondaryTrack;
67
68 for(G4int DSecLoop = 0; DSecLoop < fpParticleChange->GetNumberOfSecondaries();
69 DSecLoop++)
70 {
71 tempSecondaryTrack = fpParticleChange->GetSecondary(DSecLoop);
72
73 if(tempSecondaryTrack->GetDefinition()->GetApplyCutsFlag())
74 {
75 ApplyProductionCut(tempSecondaryTrack);
76 }
77
78 // Set parentID
79 tempSecondaryTrack->SetParentID(fpTrack->GetTrackID());
80
81 // Set the process pointer which created this track
82 tempSecondaryTrack->SetCreatorProcess(fpCurrentProcess);
83
84 // If this 2ndry particle has 'zero' kinetic energy, make sure
85 // it invokes a rest process at the beginning of the tracking
86 if(tempSecondaryTrack->GetKineticEnergy() <= DBL_MIN)
87 {
88 G4ProcessManager* pm = tempSecondaryTrack->GetDefinition()->GetProcessManager();
89 if (pm->GetAtRestProcessVector()->entries()>0)
90 {
91 tempSecondaryTrack->SetTrackStatus( fStopButAlive );
92 fpSecondary->push_back( tempSecondaryTrack );
93 fN2ndariesAtRestDoIt++;
94 }
95 else
96 {
97 delete tempSecondaryTrack;
98 }
99 }
100 else
101 {
102 fpSecondary->push_back( tempSecondaryTrack );
103 counter++;
104 }
105 } //end of loop on secondary
106 }
107
108 //_________________________________________________________________________
109
110 void G4ITStepProcessor::DoIt(double timeStep)
111
112 // Call the process having the min step length or just propagate the track on the given time step
113
114 // If the track is "leading the step" (ie one of its process has been selected
115 // as the one having the minimum time step over all tracks and processes),
116 // it will undergo its selected processes. Otherwise, it will just propagate the track
117 // on the given time step.
118
119 {
120 if(fpVerbose != nullptr) fpVerbose->DoItStarted();
121
122 G4TrackManyList* mainList = fpTrackContainer->GetMainList();
123 G4TrackManyList::iterator it = mainList->end();
124 it--;
125 std::size_t initialSize = mainList->size();
126
127 // G4cout << "initialSize = " << initialSize << G4endl;
128
129 for(std::size_t i = 0 ; i < initialSize ; ++i)
130 {
131
132 // G4cout << "i = " << i << G4endl;
133
134 G4Track* track = *it;
135 if (track == nullptr)
136 {
137 G4ExceptionDescription exceptionDescription;
138 exceptionDescription << "No track was pop back the main track list.";
139 G4Exception("G4ITStepProcessor::DoIt", "NO_TRACK",
140 FatalException, exceptionDescription);
141 }
142 // G4TrackManyList::iterator next_it (it);
143 // next_it--;
144 // it = next_it;
145
146 it--;
147 // Must be called before EndTracking(track)
148 // Otherwise the iterator will point to the list of killed tracks
149
150 if(track->GetTrackStatus() == fStopAndKill)
151 {
152 fpTrackingManager->EndTracking(track);
153 // G4cout << GetIT(track)->GetName() << G4endl;
154 // G4cout << " ************************ CONTINUE ********************" << G4endl;
155 continue;
156 }
157
158 #if defined (DEBUG_MEM) && defined (DEBUG_MEM_DOIT)
159 MemStat mem_first, mem_second, mem_diff;
160 mem_first = MemoryUsage();
161 #endif
162
163 Stepping(track, timeStep);
164
165 #if defined (DEBUG_MEM) && defined (DEBUG_MEM_DOIT)
166 MemStat mem_intermediaire = MemoryUsage();
167 mem_diff = mem_intermediaire-mem_first;
168 G4cout << "\t\t >> || MEM || In DoIT with track "
169 << track->GetTrackID() << ", diff is : " << mem_diff << G4endl;
170 #endif
171
172 ExtractDoItData();
173
174 #if defined (DEBUG_MEM) && defined (DEBUG_MEM_DOIT)
175 mem_second = MemoryUsage();
176 mem_diff = mem_second-mem_first;
177 G4cout << "\t >> || MEM || In DoIT with track "
178 << track->GetTrackID()
179 << ", diff is : " << mem_diff << G4endl;
180 #endif
181 }
182
183
184 fpTrackContainer->MergeSecondariesWithMainList();
185 fpTrackContainer->KillTracks(); // (18-06-15 : MK) Remove it ?
186 fLeadingTracks.Reset();
187 }
188
189 //_________________________________________________________________________
190
191 void G4ITStepProcessor::ExtractDoItData()
192 {
193 if (fpTrack == nullptr)
194 {
195 CleanProcessor();
196 return;
197 }
198
199 G4TrackStatus status = fpTrack->GetTrackStatus();
200
201 switch (status)
202 {
203 case fAlive:
204 case fStopButAlive:
205 case fSuspend:
206 case fPostponeToNextEvent:
207 default:
208 PushSecondaries();
209 break;
210
211 case fStopAndKill:
212 G4ITReactionSet::Instance()->RemoveReactionSet(fpTrack);
213 PushSecondaries();
214 // G4TrackList::Pop(fpTrack);
215 fpTrackingManager->EndTracking(fpTrack);
216 // fTrackContainer->PushToKill(fpTrack);
217 break;
218
219 case fKillTrackAndSecondaries:
220 G4ITReactionSet::Instance()->RemoveReactionSet(fpTrack);
221 if (fpSecondary != nullptr)
222 {
223 for (auto & i : *fpSecondary)
224 {
225 delete i;
226 }
227 fpSecondary->clear();
228 }
229 // G4TrackList::Pop(fpTrack);
230 fpTrackingManager->EndTracking(fpTrack);
231 // fTrackContainer->PushToKill(fpTrack);
232 break;
233 }
234
235 CleanProcessor();
236 }
237
238 //_________________________________________________________________________
239
240 void G4ITStepProcessor::PushSecondaries()
241 {
242 if ((fpSecondary == nullptr) || fpSecondary->empty())
243 {
244 // DEBUG
245 // G4cout << "NO SECONDARIES !!! " << G4endl;
246 return;
247 }
248
249 // DEBUG
250 // G4cout << "There are secondaries : "<< secondaries -> size() << G4endl ;
251
252 auto secondaries_i = fpSecondary->begin();
253
254 for (; secondaries_i != fpSecondary->end(); ++secondaries_i)
255 {
256 G4Track* secondary = *secondaries_i;
257 fpTrackContainer->_PushTrack(secondary);
258 }
259 }
260
261 //______________________________________________________________________________
262
263 void G4ITStepProcessor::Stepping(G4Track* track, const double & timeStep)
264 {
265
266 #ifdef DEBUG_MEM
267 MemStat mem_first, mem_second, mem_diff;
268 #endif
269
270 #ifdef DEBUG_MEM
271 mem_first = MemoryUsage();
272 #endif
273
274 CleanProcessor();
275
276 #ifdef DEBUG_MEM
277 MemStat mem_intermediaire = MemoryUsage();
278 mem_diff = mem_intermediaire-mem_first;
279 G4cout << "\t\t\t >> || MEM || After CleanProcessor " << track->GetTrackID() << ", diff is : " << mem_diff << G4endl;
280 #endif
281
282 if(track == nullptr) return; // maybe put an exception here
283 fTimeStep = timeStep;
284 SetTrack(track);
285 DoStepping();
286 }
287 //______________________________________________________________________________
288
289 // ************************************************************************
290 // Stepping
291 // ************************************************************************
292 void G4ITStepProcessor::DoStepping()
293 {
294 SetupMembers();
295
296 #ifdef DEBUG_MEM
297 MemStat mem_first, mem_second, mem_diff;
298 #endif
299
300 #ifdef DEBUG_MEM
301 mem_first = MemoryUsage();
302 #endif
303
304 #ifdef G4VERBOSE
305 if(fpVerbose != nullptr) fpVerbose->PreStepVerbose(fpTrack);
306 #endif
307
308 if(fpProcessInfo == nullptr)
309 {
310 G4ExceptionDescription exceptionDescription;
311 exceptionDescription << "No process info found for particle :"
312 << fpTrack->GetDefinition()->GetParticleName();
313 G4Exception("G4ITStepProcessor::DoStepping",
314 "ITStepProcessor0012",
315 FatalErrorInArgument,
316 exceptionDescription);
317 return;
318 }
319 // else if(fpTrack->GetTrackStatus() == fStopAndKill)
320 // {
321 // fpState->fStepStatus = fUndefined;
322 // return;
323 // }
324
325 if(fpProcessInfo->MAXofPostStepLoops == 0 &&
326 fpProcessInfo->MAXofAlongStepLoops == 0
327 && fpProcessInfo->MAXofAtRestLoops == 0)
328 {/*
329 G4ExceptionDescription exceptionDescription;
330 exceptionDescription << "No process was found for particle :"
331 << fpTrack->GetDefinition()->GetParticleName();
332 G4Exception("G4ITStepProcessor::DoStepping",
333 "ITStepProcessorNoProcess",
334 JustWarning,
335 exceptionDescription);
336
337 fpTrack->SetTrackStatus(fStopAndKill);
338 fpState->fStepStatus = fUndefined;*/
339 return;
340 }
341
342 //--------
343 // Prelude
344 //--------
345 #ifdef G4VERBOSE
346 // !!!!! Verbose
347 if(fpVerbose != nullptr) fpVerbose->NewStep();
348 #endif
349
350 //---------------------------------
351 // AtRestStep, AlongStep and PostStep Processes
352 //---------------------------------
353
354 fpNavigator->SetNavigatorState(fpITrack->GetTrackingInfo()->GetNavigatorState());
355 // fpNavigator->ResetHierarchyAndLocate( fpTrack->GetPosition(),
356 // fpTrack->GetMomentumDirection(),
357 // *((G4TouchableHistory*)fpTrack->GetTouchableHandle()()) );
358 // fpNavigator->SetNavigatorState(fpITrack->GetTrackingInfo()->GetNavigatorState());
359 // We reset the navigator state before checking for AtRest
360 // in case a AtRest processe would use a navigator info
361
362 #ifdef DEBUG_MEM
363 MemStat mem_intermediaire = MemoryUsage();
364 mem_diff = mem_intermediaire-mem_first;
365 G4cout << "\t\t\t >> || MEM || G4ITStepProcessor::DoStepping || After dealing with navigator with " << fpTrack->GetTrackID() << ", diff is : " << mem_diff << G4endl;
366 #endif
367
368 if(fpTrack->GetTrackStatus() == fStopButAlive)
369 {
370 if(fpProcessInfo->MAXofAtRestLoops > 0 && fpProcessInfo->fpAtRestDoItVector
371 != nullptr) // second condition to make coverity happy
372 {
373 //-----------------
374 // AtRestStepDoIt
375 //-----------------
376 InvokeAtRestDoItProcs();
377 fpState->fStepStatus = fAtRestDoItProc;
378 fpStep->GetPostStepPoint()->SetStepStatus(fpState->fStepStatus);
379
380 #ifdef G4VERBOSE
381 // !!!!! Verbose
382 if(fpVerbose != nullptr) fpVerbose->AtRestDoItInvoked();
383 #endif
384
385 }
386 // Make sure the track is killed
387 // fpTrack->SetTrackStatus(fStopAndKill);
388 }
389 else // if(fTimeStep > 0.) // Bye, because PostStepIL can return 0 => time =0
390 {
391 if(fpITrack == nullptr)
392 {
393 G4ExceptionDescription exceptionDescription;
394 exceptionDescription << " !!! TrackID : " << fpTrack->GetTrackID()
395 << G4endl<< " !!! Track status : "<< fpTrack->GetTrackStatus() << G4endl
396 << " !!! Particle Name : "<< fpTrack -> GetDefinition() -> GetParticleName() << G4endl
397 << "No G4ITStepProcessor::fpITrack found" << G4endl;
398
399 G4Exception("G4ITStepProcessor::DoStepping",
400 "ITStepProcessor0013",
401 FatalErrorInArgument,
402 exceptionDescription);
403 return; // to make coverity happy
404 }
405
406 if(!fpITrack->GetTrackingInfo()->IsLeadingStep())
407 {
408 // In case the track has NOT the minimum step length
409 // Given the final step time, the transportation
410 // will compute the final position of the particle
411 fpState->fStepStatus = fPostStepDoItProc;
412 fpStep->GetPostStepPoint()->SetProcessDefinedStep(fpTransportation);
413 FindTransportationStep();
414 }
415
416 #ifdef DEBUG_MEM
417 mem_intermediaire = MemoryUsage();
418 mem_diff = mem_intermediaire-mem_first;
419 G4cout << "\t\t\t >> || MEM || G4ITStepProcessor::DoStepping || After FindTransportationStep() with " << fpTrack->GetTrackID() << ", diff is : " << mem_diff << G4endl;
420 #endif
421
422 // Store the Step length (geometrical length) to G4Step and G4Track
423 fpTrack->SetStepLength(fpState->fPhysicalStep);
424 fpStep->SetStepLength(fpState->fPhysicalStep);
425
426 G4double GeomStepLength = fpState->fPhysicalStep;
427
428 // Store StepStatus to PostStepPoint
429 fpStep->GetPostStepPoint()->SetStepStatus(fpState->fStepStatus);
430
431 // Invoke AlongStepDoIt
432 InvokeAlongStepDoItProcs();
433
434 #ifdef DEBUG_MEM
435 mem_intermediaire = MemoryUsage();
436 mem_diff = mem_intermediaire-mem_first;
437 G4cout << "\t\t\t >> || MEM || G4ITStepProcessor::DoStepping || After InvokeAlongStepDoItProcs() with " << fpTrack->GetTrackID() << ", diff is : " << mem_diff << G4endl;
438 #endif
439
440 #ifdef G4VERBOSE
441 // !!!!! Verbose
442 if(fpVerbose != nullptr) fpVerbose->AlongStepDoItAllDone();
443 #endif
444
445 // Update track by taking into account all changes by AlongStepDoIt
446 // fpStep->UpdateTrack(); // done in InvokeAlongStepDoItProcs
447
448 // Update safety after invocation of all AlongStepDoIts
449 fpState->fEndpointSafOrigin = fpPostStepPoint->GetPosition();
450
451 fpState->fEndpointSafety =
452 std::max(fpState->fProposedSafety - GeomStepLength, kCarTolerance);
453
454 fpStep->GetPostStepPoint()->SetSafety(fpState->fEndpointSafety);
455
456 if(GetIT(fpTrack)->GetTrackingInfo()->IsLeadingStep())
457 {
458 // Invoke PostStepDoIt including G4ITTransportation::PSDI
459 InvokePostStepDoItProcs();
460
461 #ifdef DEBUG_MEM
462 mem_intermediaire = MemoryUsage();
463 mem_diff = mem_intermediaire-mem_first;
464 G4cout << "\t\t\t >> || MEM || G4ITStepProcessor::DoStepping || After InvokePostStepDoItProcs() with " << fpTrack->GetTrackID() << ", diff is : " << mem_diff << G4endl;
465 #endif
466 #ifdef G4VERBOSE
467 // !!!!! Verbose
468 if(fpVerbose != nullptr) fpVerbose->StepInfoForLeadingTrack();
469 #endif
470 }
471 else
472 {
473 // Only invoke transportation and all other forced processes
474 InvokeTransportationProc();
475 fpStep->GetPostStepPoint()->SetProcessDefinedStep(fpTransportation);
476
477 #ifdef DEBUG_MEM
478 mem_intermediaire = MemoryUsage();
479 mem_diff = mem_intermediaire-mem_first;
480 G4cout << "\t\t\t >> || MEM || G4ITStepProcessor::DoStepping || After InvokeTransportationProc() with " << fpTrack->GetTrackID() << ", diff is : " << mem_diff << G4endl;
481 #endif
482 }
483
484 #ifdef G4VERBOSE
485 // !!!!! Verbose
486 if(fpVerbose != nullptr) fpVerbose->PostStepDoItAllDone();
487 #endif
488 }
489
490 fpNavigator->ResetNavigatorState();
491
492 #ifdef DEBUG_MEM
493 mem_intermediaire = MemoryUsage();
494 mem_diff = mem_intermediaire-mem_first;
495 G4cout << "\t\t\t >> || MEM || G4ITStepProcessor::DoStepping || After fpNavigator->SetNavigatorState with " << fpTrack->GetTrackID() << ", diff is : " << mem_diff << G4endl;
496 #endif
497
498 //-------
499 // Finale
500 //-------
501
502 // Update 'TrackLength' and remeber the Step length of the current Step
503 fpTrack->AddTrackLength(fpStep->GetStepLength());
504 fpTrack->IncrementCurrentStepNumber();
505
506 //#ifdef G4VERBOSE
507 // // !!!!! Verbose
508 // if(fpVerbose) fpVerbose->StepInfo();
509 //#endif
510
511 #ifdef G4VERBOSE
512 if(fpVerbose != nullptr) fpVerbose->PostStepVerbose(fpTrack);
513 #endif
514
515 // G4cout << " G4ITStepProcessor::DoStepping -- " <<fpTrack->GetTrackID() << " tps = " << fpTrack->GetGlobalTime() << G4endl;
516
517 // Send G4Step information to Hit/Dig if the volume is sensitive
518 /***
519 fpCurrentVolume = fpStep->GetPreStepPoint()->GetPhysicalVolume();
520 StepControlFlag = fpStep->GetControlFlag();
521
522 if( fpCurrentVolume != 0 && StepControlFlag != AvoidHitInvocation)
523 {
524 fpSensitive = fpStep->GetPreStepPoint()->
525 GetSensitiveDetector();
526 if( fpSensitive != 0 )
527 {
528 fpSensitive->Hit(fpStep);
529 }
530 }
531
532 User intervention process.
533 if( fpUserSteppingAction != 0 )
534 {
535 fpUserSteppingAction->UserSteppingAction(fpStep);
536 }
537 G4UserSteppingAction* regionalAction
538 = fpStep->GetPreStepPoint()->GetPhysicalVolume()->GetLogicalVolume()->GetRegion()
539 ->GetRegionalSteppingAction();
540 if( regionalAction ) regionalAction->UserSteppingAction(fpStep);
541 ***/
542 fpTrackingManager->AppendStep(fpTrack, fpStep);
543 // Stepping process finish. Return the value of the StepStatus.
544
545 #ifdef DEBUG_MEM
546 MemStat mem_intermediaire = MemoryUsage();
547 mem_diff = mem_intermediaire-mem_first;
548 G4cout << "\t\t\t >> || MEM || End of DoStepping() with " << fpTrack->GetTrackID() << ", diff is : " << mem_diff << G4endl;
549 #endif
550
551 // return fpState->fStepStatus;
552 }
553
554 //______________________________________________________________________________
555
556 // ************************************************************************
557 // AtRestDoIt
558 // ************************************************************************
559
560 void G4ITStepProcessor::InvokeAtRestDoItProcs()
561 {
562 fpStep->SetStepLength(0.); //the particle has stopped
563 fpTrack->SetStepLength(0.);
564
565 G4SelectedAtRestDoItVector& selectedAtRestDoItVector =
566 fpState->fSelectedAtRestDoItVector;
567
568 // invoke selected process
569 for(std::size_t np = 0; np < fpProcessInfo->MAXofAtRestLoops; ++np)
570 {
571 //
572 // Note: DoItVector has inverse order against GetPhysIntVector
573 // and SelectedAtRestDoItVector.
574 //
575 if(selectedAtRestDoItVector[fpProcessInfo->MAXofAtRestLoops - np - 1] != InActivated)
576 {
577 fpCurrentProcess =
578 (G4VITProcess*) (*fpProcessInfo->fpAtRestDoItVector)[(G4int)np];
579
580 // G4cout << " Invoke : "
581 // << fpCurrentProcess->GetProcessName()
582 // << G4endl;
583
584 // if(fpVerbose)
585 // {
586 // fpVerbose->AtRestDoItOneByOne();
587 // }
588
589 fpCurrentProcess->SetProcessState(fpTrackingInfo->GetProcessState(fpCurrentProcess
590 ->GetProcessID()));
591 fpParticleChange = fpCurrentProcess->AtRestDoIt(*fpTrack, *fpStep);
592 fpCurrentProcess->ResetProcessState();
593
594 // Set the current process as a process which defined this Step length
595 fpStep->GetPostStepPoint()->SetProcessDefinedStep(fpCurrentProcess);
596
597 // Update Step
598 fpParticleChange->UpdateStepForAtRest(fpStep);
599
600 // Now Store the secondaries from ParticleChange to SecondaryList
601 DealWithSecondaries(fN2ndariesAtRestDoIt);
602
603 // Set the track status according to what the process defined
604 // if kinetic energy >0, otherwise set fStopButAlive
605 fpTrack->SetTrackStatus(fpParticleChange->GetTrackStatus());
606
607 // clear ParticleChange
608 fpParticleChange->Clear();
609
610 } //if(fSelectedAtRestDoItVector[np] != InActivated){
611 } //for(std::size_t np=0; np < MAXofAtRestLoops; ++np){
612 fpStep->UpdateTrack();
613
614 // Modification par rapport au transport standard :
615 // fStopAndKill doit etre propose par le modele
616 // sinon d autres processus AtRest seront appeles
617 // au pas suivant
618 // fpTrack->SetTrackStatus(fStopAndKill);
619 }
620
621 //______________________________________________________________________________
622
623 // ************************************************************************
624 // AlongStepDoIt
625 // ************************************************************************
626
627 void G4ITStepProcessor::InvokeAlongStepDoItProcs()
628 {
629
630 #ifdef DEBUG_MEM
631 MemStat mem_first, mem_second, mem_diff;
632 #endif
633
634 #ifdef DEBUG_MEM
635 mem_first = MemoryUsage();
636 #endif
637
638 // If the current Step is defined by a 'ExclusivelyForced'
639 // PostStepDoIt, then don't invoke any AlongStepDoIt
640 if(fpState->fStepStatus == fExclusivelyForcedProc)
641 {
642 return; // Take note 'return' at here !!!
643 }
644
645 // Invoke the all active continuous processes
646 for(std::size_t ci = 0; ci < fpProcessInfo->MAXofAlongStepLoops; ++ci)
647 {
648 fpCurrentProcess =
649 (G4VITProcess*) (*fpProcessInfo->fpAlongStepDoItVector)[(G4int)ci];
650 if(fpCurrentProcess == nullptr) continue;
651 // NULL means the process is inactivated by a user on fly.
652
653 fpCurrentProcess->SetProcessState(fpTrackingInfo->GetProcessState(fpCurrentProcess
654 ->GetProcessID()));
655 fpParticleChange = fpCurrentProcess->AlongStepDoIt(*fpTrack, *fpStep);
656
657 #ifdef DEBUG_MEM
658 MemStat mem_intermediaire = MemoryUsage();
659 mem_diff = mem_intermediaire-mem_first;
660 G4cout << "\t\t\t >> || MEM || After calling AlongStepDoIt for " << fpCurrentProcess->GetProcessName() << " and track "<< fpTrack->GetTrackID() << ", diff is : " << mem_diff << G4endl;
661 #endif
662
663 // fpCurrentProcess->SetProcessState(0);
664 fpCurrentProcess->ResetProcessState();
665 // Update the PostStepPoint of Step according to ParticleChange
666
667 fpParticleChange->UpdateStepForAlongStep(fpStep);
668
669 #ifdef G4VERBOSE
670 // !!!!! Verbose
671 if(fpVerbose != nullptr) fpVerbose->AlongStepDoItOneByOne();
672 #endif
673
674 // Now Store the secondaries from ParticleChange to SecondaryList
675 DealWithSecondaries(fN2ndariesAlongStepDoIt);
676
677 // Set the track status according to what the process defined
678 // if kinetic energy >0, otherwise set fStopButAlive
679 fpTrack->SetTrackStatus(fpParticleChange->GetTrackStatus());
680
681 // clear ParticleChange
682 fpParticleChange->Clear();
683 }
684
685 #ifdef DEBUG_MEM
686 MemStat mem_intermediaire = MemoryUsage();
687 mem_diff = mem_intermediaire-mem_first;
688 G4cout << "\t\t\t >> || MEM || After looping on processes with " << fpTrack->GetTrackID() << ", diff is : " << mem_diff << G4endl;
689 #endif
690
691 fpStep->UpdateTrack();
692
693 G4TrackStatus fNewStatus = fpTrack->GetTrackStatus();
694
695 if(fNewStatus == fAlive && fpTrack->GetKineticEnergy() <= DBL_MIN)
696 {
697 // G4cout << "G4ITStepProcessor::InvokeAlongStepDoItProcs : Track will be killed" << G4endl;
698 if(fpProcessInfo->MAXofAtRestLoops>0) fNewStatus = fStopButAlive;
699 else fNewStatus = fStopAndKill;
700 fpTrack->SetTrackStatus( fNewStatus );
701 }
702
703 }
704
705 //______________________________________________________________________________
706
707 // ************************************************************************
708 // PostStepDoIt
709 // ************************************************************************
710
711 void G4ITStepProcessor::InvokePostStepDoItProcs()
712 {
713 std::size_t _MAXofPostStepLoops = fpProcessInfo->MAXofPostStepLoops;
714 G4SelectedPostStepDoItVector& selectedPostStepDoItVector = fpState
715 ->fSelectedPostStepDoItVector;
716 G4StepStatus& stepStatus = fpState->fStepStatus;
717
718 // Invoke the specified discrete processes
719 for(std::size_t np = 0; np < _MAXofPostStepLoops; ++np)
720 {
721 //
722 // Note: DoItVector has inverse order against GetPhysIntVector
723 // and SelectedPostStepDoItVector.
724 //
725 G4int Cond = selectedPostStepDoItVector[_MAXofPostStepLoops
726 - np - 1];
727 if(Cond != InActivated)
728 {
729 if(((Cond == NotForced) && (stepStatus == fPostStepDoItProc)) ||
730 ((Cond == Forced) && (stepStatus != fExclusivelyForcedProc))
731 ||
732 // ((Cond == Conditionally) && (stepStatus == fAlongStepDoItProc)) ||
733 ((Cond == ExclusivelyForced) && (stepStatus == fExclusivelyForcedProc))
734 || ((Cond == StronglyForced)))
735 {
736
737 InvokePSDIP(np);
738 }
739 } //if(*fSelectedPostStepDoItVector(np)........
740
741 // Exit from PostStepLoop if the track has been killed,
742 // but extra treatment for processes with Strongly Forced flag
743 if(fpTrack->GetTrackStatus() == fStopAndKill)
744 {
745 for(std::size_t np1 = np + 1; np1 < _MAXofPostStepLoops; ++np1)
746 {
747 G4int Cond2 = selectedPostStepDoItVector[_MAXofPostStepLoops
748 - np1 - 1];
749 if(Cond2 == StronglyForced)
750 {
751 InvokePSDIP(np1);
752 }
753 }
754 break;
755 }
756 } //for(std::size_t np=0; np < MAXofPostStepLoops; ++np){
757 }
758
759 //______________________________________________________________________________
760
761 void G4ITStepProcessor::InvokePSDIP(std::size_t np)
762 {
763 fpCurrentProcess = (G4VITProcess*) (*fpProcessInfo->fpPostStepDoItVector)[(G4int)np];
764
765 fpCurrentProcess->SetProcessState(fpTrackingInfo->GetProcessState(fpCurrentProcess
766 ->GetProcessID()));
767 fpParticleChange = fpCurrentProcess->PostStepDoIt(*fpTrack, *fpStep);
768 // fpCurrentProcess->SetProcessState(0);
769 fpCurrentProcess->ResetProcessState();
770
771 // Update PostStepPoint of Step according to ParticleChange
772 fpParticleChange->UpdateStepForPostStep(fpStep);
773
774 #ifdef G4VERBOSE
775 // !!!!! Verbose
776 if(fpVerbose != nullptr) fpVerbose->PostStepDoItOneByOne();
777 #endif
778
779 // Update G4Track according to ParticleChange after each PostStepDoIt
780 fpStep->UpdateTrack();
781
782 // Update safety after each invocation of PostStepDoIts
783 fpStep->GetPostStepPoint()->SetSafety(CalculateSafety());
784
785 // Now Store the secondaries from ParticleChange to SecondaryList
786 DealWithSecondaries(fN2ndariesPostStepDoIt);
787
788 // Set the track status according to what the process defined
789 fpTrack->SetTrackStatus(fpParticleChange->GetTrackStatus());
790
791 // clear ParticleChange
792 fpParticleChange->Clear();
793 }
794
795 //______________________________________________________________________________
796
797 // ************************************************************************
798 // Transport on a given time
799 // ************************************************************************
800
801 void G4ITStepProcessor::FindTransportationStep()
802 {
803 double physicalStep(0.);
804
805 fpTransportation = fpProcessInfo->fpTransportation;
806 // dynamic_cast<G4ITTransportation*>((fpProcessInfo->fpAlongStepGetPhysIntVector)[MAXofAlongStepLoops-1]);
807
808 if(fpTrack == nullptr)
809 {
810 G4ExceptionDescription exceptionDescription;
811 exceptionDescription << "No G4ITStepProcessor::fpTrack found";
812 G4Exception("G4ITStepProcessor::FindTransportationStep",
813 "ITStepProcessor0013",
814 FatalErrorInArgument,
815 exceptionDescription);
816 return;
817
818 }
819 if(fpITrack == nullptr)
820 {
821 G4ExceptionDescription exceptionDescription;
822 exceptionDescription << "No G4ITStepProcessor::fITrack";
823 G4Exception("G4ITStepProcessor::FindTransportationStep",
824 "ITStepProcessor0014",
825 FatalErrorInArgument,
826 exceptionDescription);
827 return;
828 }
829 if((fpITrack->GetTrack()) == nullptr)
830 {
831 G4ExceptionDescription exceptionDescription;
832 exceptionDescription << "No G4ITStepProcessor::fITrack->GetTrack()";
833 G4Exception("G4ITStepProcessor::FindTransportationStep",
834 "ITStepProcessor0015",
835 FatalErrorInArgument,
836 exceptionDescription);
837 return;
838 }
839
840 if(fpTransportation != nullptr)
841 {
842 fpTransportation->SetProcessState(fpTrackingInfo->GetProcessState(fpTransportation
843 ->GetProcessID()));
844 fpTransportation->ComputeStep(*fpTrack, *fpStep, fTimeStep, physicalStep);
845
846 // fpTransportation->SetProcessState(0);
847 fpTransportation->ResetProcessState();
848 }
849
850 if(physicalStep >= DBL_MAX)
851 {
852 // G4cout << "---- 2) Setting stop and kill for " << GetIT(fpTrack)->GetName() << G4endl;
853 fpTrack -> SetTrackStatus(fStopAndKill);
854 return;
855 }
856
857 fpState->fPhysicalStep = physicalStep;
858 }
859
860 //______________________________________________________________________________
861
862 void G4ITStepProcessor::InvokeTransportationProc()
863 {
864 std::size_t _MAXofPostStepLoops = fpProcessInfo->MAXofPostStepLoops;
865 G4SelectedPostStepDoItVector& selectedPostStepDoItVector = fpState
866 ->fSelectedPostStepDoItVector;
867 G4StepStatus& stepStatus = fpState->fStepStatus;
868
869 // Invoke the specified discrete processes
870 for(std::size_t np = 0; np < _MAXofPostStepLoops; ++np)
871 {
872 //
873 // Note: DoItVector has inverse order against GetPhysIntVector
874 // and SelectedPostStepDoItVector.
875 //
876 G4int Cond = selectedPostStepDoItVector[_MAXofPostStepLoops - np - 1];
877 if(Cond != InActivated)
878 {
879 if(((Cond == Forced) && (stepStatus != fExclusivelyForcedProc)) ||
880 // ((Cond == Conditionally) && (stepStatus == fAlongStepDoItProc)) ||
881 ((Cond == ExclusivelyForced) && (stepStatus == fExclusivelyForcedProc))
882 || ((Cond == StronglyForced)))
883 {
884
885 InvokePSDIP(np);
886 }
887 } //if(Cond != InActivated)
888
889 // Exit from PostStepLoop if the track has been killed,
890 // but extra treatment for processes with Strongly Forced flag
891 if(fpTrack->GetTrackStatus() == fStopAndKill)
892 {
893 for(std::size_t np1 = np + 1; np1 < _MAXofPostStepLoops; ++np1)
894 {
895 G4int Cond2 = selectedPostStepDoItVector[_MAXofPostStepLoops - np1 - 1];
896 if(Cond2 == StronglyForced)
897 {
898 InvokePSDIP(np1);
899 }
900 }
901 break;
902 }
903 }
904 }
905
906 //______________________________________________________________________________
907
908 // ************************************************************************
909 // Apply cuts
910 // ************************************************************************
911
912 void G4ITStepProcessor::ApplyProductionCut(G4Track* aSecondary)
913 {
914 G4bool tBelowCutEnergyAndSafety = false;
915 G4int tPtclIdx = G4ProductionCuts::GetIndex(aSecondary->GetDefinition());
916 if(tPtclIdx < 0)
917 {
918 return;
919 }
920 G4ProductionCutsTable* tCutsTbl =
921 G4ProductionCutsTable::GetProductionCutsTable();
922 G4int tCoupleIdx = tCutsTbl->GetCoupleIndex(fpPreStepPoint
923 ->GetMaterialCutsCouple());
924 G4double tProdThreshold =
925 (*(tCutsTbl->GetEnergyCutsVector(tPtclIdx)))[tCoupleIdx];
926 if(aSecondary->GetKineticEnergy() < tProdThreshold)
927 {
928 tBelowCutEnergyAndSafety = true;
929 if(std::abs(aSecondary->GetDynamicParticle()->GetCharge()) > DBL_MIN)
930 {
931 G4double currentRange
932 = G4LossTableManager::Instance()->GetRange(aSecondary->GetDefinition(),
933 aSecondary->GetKineticEnergy(),
934 fpPreStepPoint->GetMaterialCutsCouple());
935 tBelowCutEnergyAndSafety = (currentRange < CalculateSafety() );
936 }
937 }
938
939 if(tBelowCutEnergyAndSafety)
940 {
941 if(!(aSecondary->IsGoodForTracking()))
942 {
943 // Add kinetic energy to the total energy deposit
944 fpStep->AddTotalEnergyDeposit(aSecondary->GetKineticEnergy());
945 aSecondary->SetKineticEnergy(0.0);
946 }
947 }
948 }
949