| 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
37 #include "G4ForceCondition.hh"
38 #include "G4GPILSelection.hh"
39 #include "G4GeometryTolerance.hh"
40 #include "G4IT.hh"
41 #include "G4ITNavigator.hh" // Include from 'geometry'
42 #include "G4ITSteppingVerbose.hh"
43 #include "G4ITTrackHolder.hh"
44 #include "G4ITTrackingInteractivity.hh"
45 #include "G4ITTrackingManager.hh"
46 #include "G4ITTransportation.hh"
47 #include "G4ITTransportationManager.hh"
48 #include "G4ParticleTable.hh"
49 #include "G4TrackingInformation.hh"
50 #include "G4UImanager.hh"
51 #include "G4VITProcess.hh"
52 #include "G4VITSteppingVerbose.hh"
53 #include "G4VProcess.hh"
54
55 #include <iomanip> // Include from 'system'
56 #include <vector> // Include from 'system'
57
58 using namespace std;
59
60 static const std::size_t SizeOfSelectedDoItVector = 100;
61
62 template<typename T>
63 inline bool IsInf(T value)
64 {
65 return std::numeric_limits<T>::has_infinity
66 && value == std::numeric_limits<T>::infinity();
67 }
68
69 //______________________________________________________________________________
70
71 G4ITStepProcessor::G4ITStepProcessor()
72 {
73 fpVerbose = nullptr;
74 // fpUserSteppingAction = 0 ;
75 fStoreTrajectory = 0;
76 fpTrackingManager = nullptr;
77 fpNavigator = nullptr;
78 kCarTolerance = -1.;
79 fInitialized = false;
80 fPreviousTimeStep = DBL_MAX;
81 fILTimeStep = DBL_MAX;
82 fpTrackContainer = nullptr;
83
84 CleanProcessor();
85 ResetSecondaries();
86 }
87
88 //______________________________________________________________________________
89
90 //G4ITStepProcessor::
91 G4ITStepProcessorState::G4ITStepProcessorState() :
92 G4ITStepProcessorState_Lock(),
93 fSelectedAtRestDoItVector(G4VITProcess::GetMaxProcessIndex(), 0),
94 fSelectedPostStepDoItVector(G4VITProcess::GetMaxProcessIndex(), 0)
95 {
96 fPhysicalStep = -1.;
97 fPreviousStepSize = -1.;
98
99 fSafety = -1.;
100 fProposedSafety = -1.;
101 fEndpointSafety = -1;
102
103 fStepStatus = fUndefined;
104
105 fTouchableHandle = nullptr;
106 }
107
108 //______________________________________________________________________________
109
110 //G4ITStepProcessor::
111 G4ITStepProcessorState::G4ITStepProcessorState(const G4ITStepProcessorState& right) :
112 G4ITStepProcessorState_Lock(),
113 fSelectedAtRestDoItVector(right.fSelectedAtRestDoItVector),
114 fSelectedPostStepDoItVector(right.fSelectedPostStepDoItVector)
115 {
116 fPhysicalStep = right.fPhysicalStep;
117 fPreviousStepSize = right.fPreviousStepSize;
118
119 fSafety = right.fSafety;
120 fProposedSafety = right.fProposedSafety;
121 fEndpointSafety = right.fEndpointSafety;
122
123 fStepStatus = right.fStepStatus;
124
125 fTouchableHandle = right.fTouchableHandle;
126 }
127
128 //______________________________________________________________________________
129
130 //G4ITStepProcessor::
131 G4ITStepProcessorState&
132 //G4ITStepProcessor::
133 G4ITStepProcessorState::operator=(const G4ITStepProcessorState& right)
134 {
135 if(this == &right) return *this;
136
137 fSelectedAtRestDoItVector.clear();
138 fSelectedAtRestDoItVector = right.fSelectedAtRestDoItVector;
139 fSelectedPostStepDoItVector.clear();
140 fSelectedPostStepDoItVector = right.fSelectedPostStepDoItVector;
141
142 fPhysicalStep = right.fPhysicalStep;
143 fPreviousStepSize = right.fPreviousStepSize;
144
145 fSafety = right.fSafety;
146 fProposedSafety = right.fProposedSafety;
147 fEndpointSafety = right.fEndpointSafety;
148
149 fStepStatus = right.fStepStatus;
150
151 fTouchableHandle = right.fTouchableHandle;
152 return *this;
153 }
154
155 //______________________________________________________________________________
156
157 //G4ITStepProcessor::
158 G4ITStepProcessorState::~G4ITStepProcessorState()
159 {
160 ;
161 }
162
163 //______________________________________________________________________________
164
165 void G4ITStepProcessor::ClearProcessInfo()
166 {
167 std::map<const G4ParticleDefinition*, ProcessGeneralInfo*>::iterator it;
168
169 for(it = fProcessGeneralInfoMap.begin(); it != fProcessGeneralInfoMap.end();
170 it++)
171 {
172 if(it->second != nullptr)
173 {
174 delete it->second;
175 it->second = 0;
176 }
177 }
178
179 fProcessGeneralInfoMap.clear();
180 }
181
182 //______________________________________________________________________________
183
184 void G4ITStepProcessor::ForceReInitialization()
185 {
186 fInitialized = false;
187 ClearProcessInfo();
188 Initialize();
189 }
190
191 //______________________________________________________________________________
192
193 void G4ITStepProcessor::Initialize()
194 {
195 CleanProcessor();
196 if(fInitialized) return;
197 // ActiveOnlyITProcess();
198
199 SetNavigator(G4ITTransportationManager::GetTransportationManager()
200 ->GetNavigatorForTracking());
201
202 fPhysIntLength = DBL_MAX;
203 kCarTolerance = 0.5
204 * G4GeometryTolerance::GetInstance()->GetSurfaceTolerance();
205
206 if(fpVerbose == nullptr)
207 {
208 G4ITTrackingInteractivity* interactivity = fpTrackingManager->GetInteractivity();
209
210 if(interactivity != nullptr)
211 {
212 fpVerbose = interactivity->GetSteppingVerbose();
213 fpVerbose->SetStepProcessor(this);
214 }
215 }
216
217 fpTrackContainer = G4ITTrackHolder::Instance();
218
219 fInitialized = true;
220 }
221 //______________________________________________________________________________
222
223 G4ITStepProcessor::~G4ITStepProcessor()
224 {
225 if(fpStep != nullptr)
226 {
227 fpStep->DeleteSecondaryVector();
228 delete fpStep;
229 }
230
231 delete fpSecondary;
232 ClearProcessInfo();
233 //G4ITTransportationManager::DeleteInstance();
234
235 // if(fpUserSteppingAction) delete fpUserSteppingAction;
236 }
237 //______________________________________________________________________________
238 // should not be used
239 G4ITStepProcessor::G4ITStepProcessor(const G4ITStepProcessor& rhs)
240 {
241 fpVerbose = rhs.fpVerbose;
242 fStoreTrajectory = rhs.fStoreTrajectory;
243
244 // fpUserSteppingAction = 0 ;
245 fpTrackingManager = nullptr;
246 fpNavigator = nullptr;
247 fInitialized = false;
248
249 kCarTolerance = rhs.kCarTolerance;
250 fInitialized = false;
251 fPreviousTimeStep = DBL_MAX;
252
253 CleanProcessor();
254 ResetSecondaries();
255 fpTrackContainer = nullptr;
256 fILTimeStep = DBL_MAX;
257 }
258
259 //______________________________________________________________________________
260
261 void G4ITStepProcessor::ResetLeadingTracks()
262 {
263 fLeadingTracks.Reset();
264 }
265
266 //______________________________________________________________________________
267
268 void G4ITStepProcessor::PrepareLeadingTracks()
269 {
270 fLeadingTracks.PrepareLeadingTracks();
271 }
272
273 //______________________________________________________________________________
274
275 G4ITStepProcessor& G4ITStepProcessor::operator=(const G4ITStepProcessor& rhs)
276 {
277 if(this == &rhs) return *this; // handle self assignment
278 //assignment operator
279 return *this;
280 }
281
282 //______________________________________________________________________________
283
284 void G4ITStepProcessor::ActiveOnlyITProcess()
285 {
286 // Method not used for the time being
287 #ifdef debug
288 G4cout<<"G4ITStepProcessor::CloneProcesses: is called"<<G4endl;
289 #endif
290
291 G4ParticleTable* theParticleTable = G4ParticleTable::GetParticleTable();
292 G4ParticleTable::G4PTblDicIterator* theParticleIterator = theParticleTable
293 ->GetIterator();
294
295 theParticleIterator->reset();
296 // TODO : Ne faire la boucle que sur les IT **** !!!
297 while((*theParticleIterator)())
298 {
299 G4ParticleDefinition* particle = theParticleIterator->value();
300 G4ProcessManager* pm = particle->GetProcessManager();
301
302 if(pm == nullptr)
303 {
304 G4cerr << "ERROR - G4ITStepProcessor::GetProcessNumber()" << G4endl<< " ProcessManager is NULL for particle = "
305 << particle->GetParticleName() << ", PDG_code = "
306 << particle->GetPDGEncoding() << G4endl;
307 G4Exception("G4ITStepProcessor::GetProcessNumber()", "ITStepProcessor0001",
308 FatalException, "Process Manager is not found.");
309 return;
310 }
311
312 ActiveOnlyITProcess(pm);
313 }
314 }
315
316 //______________________________________________________________________________
317
318 void G4ITStepProcessor::ActiveOnlyITProcess(G4ProcessManager* processManager)
319 {
320 // Method not used for the time being
321 G4ProcessVector* processVector = processManager->GetProcessList();
322
323 G4VITProcess* itProcess = nullptr;
324 for(G4int i = 0; i < (G4int)processVector->size(); ++i)
325 {
326 G4VProcess* base_process = (*processVector)[i];
327 itProcess = dynamic_cast<G4VITProcess*>(base_process);
328
329 if(itProcess == nullptr)
330 {
331 processManager->SetProcessActivation(base_process, false);
332 }
333 }
334 }
335
336 //______________________________________________________________________________
337
338 void G4ITStepProcessor::SetupGeneralProcessInfo(G4ParticleDefinition* particle,
339 G4ProcessManager* pm)
340 {
341
342 #ifdef debug
343 G4cout<<"G4ITStepProcessor::GetProcessNumber: is called track"<<G4endl;
344 #endif
345 if(pm == nullptr)
346 {
347 G4cerr << "ERROR - G4SteppingManager::GetProcessNumber()" << G4endl<< " ProcessManager is NULL for particle = "
348 << particle->GetParticleName() << ", PDG_code = "
349 << particle->GetPDGEncoding() << G4endl;
350 G4Exception("G4SteppingManager::GetProcessNumber()", "ITStepProcessor0002",
351 FatalException, "Process Manager is not found.");
352 return;
353 }
354
355 auto it =
356 fProcessGeneralInfoMap.find(particle);
357 if(it != fProcessGeneralInfoMap.end())
358 {
359 G4Exception("G4SteppingManager::SetupGeneralProcessInfo()",
360 "ITStepProcessor0003",
361 FatalException, "Process info already registered.");
362 return;
363 }
364
365 // here used as temporary
366 fpProcessInfo = new ProcessGeneralInfo();
367
368 // AtRestDoits
369 fpProcessInfo->MAXofAtRestLoops = pm->GetAtRestProcessVector()->entries();
370 fpProcessInfo->fpAtRestDoItVector = pm->GetAtRestProcessVector(typeDoIt);
371 fpProcessInfo->fpAtRestGetPhysIntVector =
372 pm->GetAtRestProcessVector(typeGPIL);
373 #ifdef debug
374 G4cout << "G4ITStepProcessor::GetProcessNumber: #ofAtRest="
375 << fpProcessInfo->MAXofAtRestLoops << G4endl;
376 #endif
377
378 // AlongStepDoits
379 fpProcessInfo->MAXofAlongStepLoops =
380 pm->GetAlongStepProcessVector()->entries();
381 fpProcessInfo->fpAlongStepDoItVector =
382 pm->GetAlongStepProcessVector(typeDoIt);
383 fpProcessInfo->fpAlongStepGetPhysIntVector =
384 pm->GetAlongStepProcessVector(typeGPIL);
385 #ifdef debug
386 G4cout << "G4ITStepProcessor::GetProcessNumber:#ofAlongStp="
387 << fpProcessInfo->MAXofAlongStepLoops << G4endl;
388 #endif
389
390 // PostStepDoits
391 fpProcessInfo->MAXofPostStepLoops =
392 pm->GetPostStepProcessVector()->entries();
393 fpProcessInfo->fpPostStepDoItVector = pm->GetPostStepProcessVector(typeDoIt);
394 fpProcessInfo->fpPostStepGetPhysIntVector =
395 pm->GetPostStepProcessVector(typeGPIL);
396 #ifdef debug
397 G4cout << "G4ITStepProcessor::GetProcessNumber: #ofPostStep="
398 << fpProcessInfo->MAXofPostStepLoops << G4endl;
399 #endif
400
401 if (SizeOfSelectedDoItVector<fpProcessInfo->MAXofAtRestLoops ||
402 SizeOfSelectedDoItVector<fpProcessInfo->MAXofAlongStepLoops ||
403 SizeOfSelectedDoItVector<fpProcessInfo->MAXofPostStepLoops )
404 {
405 G4cerr << "ERROR - G4ITStepProcessor::GetProcessNumber()" << G4endl
406 << " SizeOfSelectedDoItVector= " << SizeOfSelectedDoItVector
407 << " ; is smaller then one of MAXofAtRestLoops= "
408 << fpProcessInfo->MAXofAtRestLoops << G4endl
409 << " or MAXofAlongStepLoops= " << fpProcessInfo->MAXofAlongStepLoops
410 << " or MAXofPostStepLoops= " << fpProcessInfo->MAXofPostStepLoops << G4endl;
411 G4Exception("G4ITStepProcessor::GetProcessNumber()",
412 "ITStepProcessor0004", FatalException,
413 "The array size is smaller than the actual No of processes.");
414 }
415
416 if((fpProcessInfo->fpAtRestDoItVector == nullptr) &&
417 (fpProcessInfo->fpAlongStepDoItVector == nullptr) &&
418 (fpProcessInfo->fpPostStepDoItVector == nullptr))
419 {
420 G4ExceptionDescription exceptionDescription;
421 exceptionDescription << "No DoIt process found ";
422 G4Exception("G4ITStepProcessor::DoStepping","ITStepProcessor0005",
423 FatalErrorInArgument,exceptionDescription);
424 return;
425 }
426
427 if((fpProcessInfo->fpAlongStepGetPhysIntVector != nullptr)
428 && fpProcessInfo->MAXofAlongStepLoops>0)
429 {
430 fpProcessInfo->fpTransportation = dynamic_cast<G4ITTransportation*>
431 ((*fpProcessInfo->fpAlongStepGetPhysIntVector)
432 [G4int(fpProcessInfo->MAXofAlongStepLoops-1)]);
433
434 if(fpProcessInfo->fpTransportation == nullptr)
435 {
436 G4ExceptionDescription exceptionDescription;
437 exceptionDescription << "No transportation process found ";
438 G4Exception("G4ITStepProcessor::SetupGeneralProcessInfo",
439 "ITStepProcessor0006",
440 FatalErrorInArgument,exceptionDescription);
441 }
442 }
443 fProcessGeneralInfoMap[particle] = fpProcessInfo;
444 // fpProcessInfo = 0;
445 }
446
447 //______________________________________________________________________________
448
449 void G4ITStepProcessor::SetTrack(G4Track* track)
450 {
451 fpTrack = track;
452 if(fpTrack != nullptr)
453 {
454 fpITrack = GetIT(fpTrack);
455 fpStep = const_cast<G4Step*>(fpTrack->GetStep());
456
457 if(fpITrack != nullptr)
458 {
459 fpTrackingInfo = fpITrack->GetTrackingInfo();
460 }
461 else
462 {
463 fpTrackingInfo = nullptr;
464 G4cerr << "Track ID : " << fpTrack->GetTrackID() << G4endl;
465
466 G4ExceptionDescription errMsg;
467 errMsg << "No IT pointer was attached to the track you try to process.";
468 G4Exception("G4ITStepProcessor::SetTrack",
469 "ITStepProcessor0007",
470 FatalErrorInArgument,
471 errMsg);
472 }
473 }
474 else
475 {
476 fpITrack = nullptr;
477 fpStep = nullptr;
478 }
479 }
480 //______________________________________________________________________________
481
482 void G4ITStepProcessor::GetProcessInfo()
483 {
484 G4ParticleDefinition* particle = fpTrack->GetDefinition();
485 auto it =
486 fProcessGeneralInfoMap.find(particle);
487
488 if(it == fProcessGeneralInfoMap.end())
489 {
490 SetupGeneralProcessInfo(particle,
491 fpTrack->GetDefinition()->GetProcessManager());
492 if(fpProcessInfo == nullptr)
493 {
494 G4ExceptionDescription exceptionDescription("...");
495 G4Exception("G4ITStepProcessor::GetProcessNumber",
496 "ITStepProcessor0008",
497 FatalErrorInArgument,
498 exceptionDescription);
499 return;
500 }
501 }
502 else
503 {
504 fpProcessInfo = it->second;
505 }
506 }
507
508 //______________________________________________________________________________
509
510 void G4ITStepProcessor::SetupMembers()
511 {
512 fpSecondary = fpStep->GetfSecondary();
513 fpPreStepPoint = fpStep->GetPreStepPoint();
514 fpPostStepPoint = fpStep->GetPostStepPoint();
515
516 fpState = (G4ITStepProcessorState*) fpITrack->GetTrackingInfo()
517 ->GetStepProcessorState();
518
519 GetProcessInfo();
520 ResetSecondaries();
521 }
522
523 //______________________________________________________________________________
524
525 void G4ITStepProcessor::ResetSecondaries()
526 {
527 // Reset the secondary particles
528 fN2ndariesAtRestDoIt = 0;
529 fN2ndariesAlongStepDoIt = 0;
530 fN2ndariesPostStepDoIt = 0;
531 }
532
533 //______________________________________________________________________________
534
535 void G4ITStepProcessor::GetAtRestIL()
536 {
537 // Select the rest process which has the shortest time before
538 // it is invoked. In rest processes, GPIL()
539 // returns the time before a process occurs.
540 G4double lifeTime(DBL_MAX), shortestLifeTime (DBL_MAX);
541
542 fAtRestDoItProcTriggered = 0;
543 shortestLifeTime = DBL_MAX;
544
545 unsigned int NofInactiveProc=0;
546
547 for( G4int ri=0; ri < (G4int)fpProcessInfo->MAXofAtRestLoops; ++ri )
548 {
549 fpCurrentProcess = dynamic_cast<G4VITProcess*>((*fpProcessInfo->fpAtRestGetPhysIntVector)[ri]);
550 if (fpCurrentProcess== nullptr)
551 {
552 (fpState->fSelectedAtRestDoItVector)[ri] = InActivated;
553 NofInactiveProc++;
554 continue;
555 } // NULL means the process is inactivated by a user on fly.
556
557 fCondition=NotForced;
558 fpCurrentProcess->SetProcessState(
559 fpTrackingInfo->GetProcessState(fpCurrentProcess->GetProcessID()));
560
561 lifeTime = fpCurrentProcess->AtRestGPIL( *fpTrack, &fCondition );
562 fpCurrentProcess->ResetProcessState();
563
564 if(fCondition==Forced)
565 {
566 (fpState->fSelectedAtRestDoItVector)[ri] = Forced;
567 }
568 else
569 {
570 (fpState->fSelectedAtRestDoItVector)[ri] = InActivated;
571 if(lifeTime < shortestLifeTime )
572 {
573 shortestLifeTime = lifeTime;
574 fAtRestDoItProcTriggered = G4int(ri);
575 }
576 }
577 }
578
579 (fpState->fSelectedAtRestDoItVector)[fAtRestDoItProcTriggered] = NotForced;
580
581 // G4cout << " --> Selected at rest process : "
582 // << (*fpProcessInfo->fpAtRestGetPhysIntVector)[fAtRestDoItProcTriggered]
583 // ->GetProcessName()
584 // << G4endl;
585
586 fTimeStep = shortestLifeTime;
587
588 // at least one process is necessary to destroy the particle
589 // exit with warning
590 if(NofInactiveProc==fpProcessInfo->MAXofAtRestLoops)
591 {
592 G4cerr << "ERROR - G4ITStepProcessor::InvokeAtRestDoItProcs()" << G4endl
593 << " No AtRestDoIt process is active!" << G4endl;
594 }
595 }
596
597 //_________________________________________________________________________
598 G4double G4ITStepProcessor::ComputeInteractionLength(double previousTimeStep)
599 {
600 G4TrackManyList* mainList = fpTrackContainer->GetMainList();
601 G4TrackManyList::iterator it = mainList ->begin();
602 G4TrackManyList::iterator end = mainList ->end();
603
604 SetPreviousStepTime(previousTimeStep);
605
606 fILTimeStep = DBL_MAX;
607
608 for (; it != end; )
609 {
610 G4Track * track = *it;
611
612 #ifdef DEBUG
613 G4cout << "*CIL* " << GetIT(track)->GetName()
614 << " ID: " << track->GetTrackID()
615 << " at time : " << track->GetGlobalTime()
616 << G4endl;
617 #endif
618
619 ++it;
620 DefinePhysicalStepLength(track);
621
622 ExtractILData();
623 }
624
625 return fILTimeStep;
626 }
627
628 //_________________________________________________________________________
629
630 void G4ITStepProcessor::ExtractILData()
631 {
632 assert(fpTrack != 0);
633 if (fpTrack == nullptr)
634 {
635 CleanProcessor();
636 return;
637 }
638
639 // assert(fpTrack->GetTrackStatus() != fStopAndKill);
640
641 if (fpTrack->GetTrackStatus() == fStopAndKill)
642 {
643 // trackContainer->GetMainList()->pop(fpTrack);
644 fpTrackingManager->EndTracking(fpTrack);
645 CleanProcessor();
646 return;
647 }
648
649 if (IsInf(fTimeStep))
650 {
651 // G4cout << "!!!!!!!!!!!! IS INF " << track->GetTrackID() << G4endl;
652 CleanProcessor();
653 return;
654 }
655 if (fTimeStep < fILTimeStep - DBL_EPSILON)
656 {
657 // G4cout << "!!!!!!!!!!!! TEMPS DIFFERENTS "
658 // << track->GetTrackID() << G4endl;
659
660 fLeadingTracks.Reset();
661
662 fILTimeStep = GetInteractionTime();
663
664 // G4cout << "Will set leading step to true for time :"
665 // << SP -> GetInteractionTime() << " against fTimeStep : "
666 // << G4BestUnit(fILTimeStep, "Time") << " the trackID is : " << track->GetTrackID()
667 // << G4endl;
668
669 // GetIT(fpTrack)->GetTrackingInfo()->SetLeadingStep(true);
670 fLeadingTracks.Push(fpTrack);
671 }
672 else if(fabs(fILTimeStep - fTimeStep) < DBL_EPSILON )
673 {
674
675 // G4cout << "!!!!!!!!!!!! MEME TEMPS " << track->GetTrackID() << G4endl;
676 // G4cout << "Will set leading step to true for time :"
677 // << SP -> GetInteractionTime() << " against fTimeStep : "
678 // << fTimeStep << " the trackID is : " << track->GetTrackID()<< G4endl;
679 // GetIT(fpTrack)->GetTrackingInfo()->SetLeadingStep(true);
680 fLeadingTracks.Push(fpTrack);
681 }
682 // else
683 // {
684 // G4cout << "!!!! Bigger time : " << "currentTime : "<<fILTimeStep
685 // << " proposedTime : " << SP -> GetInteractionTime() << G4endl;
686 // }
687
688 CleanProcessor();
689 }
690
691 //___________________________________________________________________________
692
693 void G4ITStepProcessor::DefinePhysicalStepLength(G4Track* track)
694 {
695 SetTrack(track);
696 DoDefinePhysicalStepLength();
697 }
698
699 //______________________________________________________________________________
700
701 void G4ITStepProcessor::SetInitialStep()
702 {
703 // DEBUG
704 // G4cout << "SetInitialStep for : " << fpITrack-> GetName() << G4endl;
705 //________________________________________________________
706 // Initialize geometry
707
708 if(!fpTrack->GetTouchableHandle())
709 {
710 //==========================================================================
711 // Create navigator state and Locate particle in geometry
712 //==========================================================================
713 /*
714 fpNavigator->NewNavigatorStateAndLocate(fpTrack->GetPosition(),
715 fpTrack->GetMomentumDirection());
716
717 fpITrack->GetTrackingInfo()->
718 SetNavigatorState(fpNavigator->GetNavigatorState());
719 */
720 fpNavigator->NewNavigatorState();
721 fpITrack->GetTrackingInfo()->SetNavigatorState(fpNavigator
722 ->GetNavigatorState());
723
724 G4ThreeVector direction = fpTrack->GetMomentumDirection();
725 fpNavigator->LocateGlobalPointAndSetup(fpTrack->GetPosition(),
726 &direction,
727 false,
728 false); // was false, false
729
730 fpState->fTouchableHandle = fpNavigator->CreateTouchableHistory();
731
732 fpTrack->SetTouchableHandle(fpState->fTouchableHandle);
733 fpTrack->SetNextTouchableHandle(fpState->fTouchableHandle);
734 }
735 else
736 {
737 fpState->fTouchableHandle = fpTrack->GetTouchableHandle();
738 fpTrack->SetNextTouchableHandle(fpState->fTouchableHandle);
739
740 //==========================================================================
741 // Create OR set navigator state
742 //==========================================================================
743
744 if(fpITrack->GetTrackingInfo()->GetNavigatorState() != nullptr)
745 {
746 fpNavigator->SetNavigatorState(fpITrack->GetTrackingInfo()
747 ->GetNavigatorState());
748 fpITrack->GetTrackingInfo()->SetNavigatorState(fpNavigator
749 ->GetNavigatorState());
750 }
751 else
752 {
753 fpNavigator->NewNavigatorState(*((G4TouchableHistory*) fpState
754 ->fTouchableHandle()));
755 fpITrack->GetTrackingInfo()->SetNavigatorState(fpNavigator
756 ->GetNavigatorState());
757 }
758
759 G4VPhysicalVolume* oldTopVolume =
760 fpTrack->GetTouchableHandle()->GetVolume();
761
762 //==========================================================================
763 // Locate particle in geometry
764 //==========================================================================
765
766 // G4VPhysicalVolume* newTopVolume =
767 // fpNavigator->LocateGlobalPointAndSetup(
768 // fpTrack->GetPosition(),
769 // &fpTrack->GetMomentumDirection(),
770 // true, false);
771
772 G4VPhysicalVolume* newTopVolume =
773 fpNavigator->ResetHierarchyAndLocate(fpTrack->GetPosition(),
774 fpTrack->GetMomentumDirection(),
775 *((G4TouchableHistory*) fpTrack
776 ->GetTouchableHandle()()));
777
778 if(newTopVolume != oldTopVolume || oldTopVolume->GetRegularStructureId()
779 == 1)
780 {
781 fpState->fTouchableHandle = fpNavigator->CreateTouchableHistory();
782 fpTrack->SetTouchableHandle(fpState->fTouchableHandle);
783 fpTrack->SetNextTouchableHandle(fpState->fTouchableHandle);
784 }
785 }
786
787 fpCurrentVolume = fpState->fTouchableHandle->GetVolume();
788
789 //________________________________________________________
790 // If the primary track has 'Suspend' or 'PostponeToNextEvent' state,
791 // set the track state to 'Alive'.
792 if((fpTrack->GetTrackStatus() == fSuspend) || (fpTrack->GetTrackStatus()
793 == fPostponeToNextEvent))
794 {
795 fpTrack->SetTrackStatus(fAlive);
796 }
797
798 //HoangTRAN: it's better to check the status here
799 if(fpTrack->GetTrackStatus() == fStopAndKill) return;
800
801 // If the primary track has 'zero' kinetic energy, set the track
802 // state to 'StopButAlive'.
803 if(fpTrack->GetKineticEnergy() <= 0.0)
804 {
805 fpTrack->SetTrackStatus(fStopButAlive);
806 }
807 //________________________________________________________
808 // Set vertex information of G4Track at here
809 if(fpTrack->GetCurrentStepNumber() == 0)
810 {
811 fpTrack->SetVertexPosition(fpTrack->GetPosition());
812 fpTrack->SetVertexMomentumDirection(fpTrack->GetMomentumDirection());
813 fpTrack->SetVertexKineticEnergy(fpTrack->GetKineticEnergy());
814 fpTrack->SetLogicalVolumeAtVertex(fpTrack->GetVolume()->GetLogicalVolume());
815 }
816 //________________________________________________________
817 // If track is already outside the world boundary, kill it
818 if(fpCurrentVolume == nullptr)
819 {
820 // If the track is a primary, stop processing
821 if(fpTrack->GetParentID() == 0)
822 {
823 G4cerr << "ERROR - G4ITStepProcessor::SetInitialStep()" << G4endl<< " Primary particle starting at - "
824 << fpTrack->GetPosition()
825 << " - is outside of the world volume." << G4endl;
826 G4Exception("G4ITStepProcessor::SetInitialStep()", "ITStepProcessor0011",
827 FatalException, "Primary vertex outside of the world!");
828 }
829
830 fpTrack->SetTrackStatus( fStopAndKill );
831 G4cout << "WARNING - G4ITStepProcessor::SetInitialStep()" << G4endl
832 << " Initial track position is outside world! - "
833 << fpTrack->GetPosition() << G4endl;
834 }
835 else
836 {
837 // Initial set up for attribues of 'Step'
838 fpStep->InitializeStep( fpTrack );
839 }
840
841 fpState->fStepStatus = fUndefined;
842 }
843 //______________________________________________________________________________
844
845 void G4ITStepProcessor::InitDefineStep()
846 {
847
848 if(fpStep == nullptr)
849 {
850 // Create new Step and give it to the track
851 fpStep = new G4Step();
852 fpTrack->SetStep(fpStep);
853 fpSecondary = fpStep->NewSecondaryVector();
854
855 // Create new state and set it in the trackingInfo
856 fpState = new G4ITStepProcessorState();
857 fpITrack->GetTrackingInfo()->SetStepProcessorState((G4ITStepProcessorState_Lock*) fpState);
858
859 SetupMembers();
860 SetInitialStep();
861
862 fpTrackingManager->StartTracking(fpTrack);
863 }
864 else
865 {
866 SetupMembers();
867
868 fpState->fPreviousStepSize = fpTrack->GetStepLength();
869 /***
870 // Send G4Step information to Hit/Dig if the volume is sensitive
871 fpCurrentVolume = fpStep->GetPreStepPoint()->GetPhysicalVolume();
872 StepControlFlag = fpStep->GetControlFlag();
873 if( fpCurrentVolume != 0 && StepControlFlag != AvoidHitInvocation)
874 {
875 fpSensitive = fpStep->GetPreStepPoint()->
876 GetSensitiveDetector();
877
878 // if( fSensitive != 0 ) {
879 // fSensitive->Hit(fStep);
880 // }
881 }
882 ***/
883 // Store last PostStepPoint to PreStepPoint, and swap current and next
884 // volume information of G4Track. Reset total energy deposit in one Step.
885 fpStep->CopyPostToPreStepPoint();
886 fpStep->ResetTotalEnergyDeposit();
887
888 //JA Set the volume before it is used (in DefineStepLength() for User Limit)
889 fpCurrentVolume = fpStep->GetPreStepPoint()->GetPhysicalVolume();
890 /*
891 G4cout << G4endl;
892 G4cout << "!!!!!!!!!!!!!!!!!!!!!!!!!!!" << G4endl;
893 G4cout << "PreStepPoint Volume : "
894 << fpCurrentVolume->GetName() << G4endl;
895 G4cout << "Track Touchable : "
896 << fpTrack->GetTouchableHandle()->GetVolume()->GetName() << G4endl;
897 G4cout << "Track NextTouchable : "
898 << fpTrack->GetNextTouchableHandle()->GetVolume()->GetName()
899 << G4endl;
900 */
901 // Reset the step's auxiliary points vector pointer
902 fpStep->SetPointerToVectorOfAuxiliaryPoints(nullptr);
903
904 // Switch next touchable in track to current one
905 fpTrack->SetTouchableHandle(fpTrack->GetNextTouchableHandle());
906 fpState->fTouchableHandle = fpTrack->GetTouchableHandle();
907 fpTrack->SetNextTouchableHandle(fpState->fTouchableHandle);
908
909 //! ADDED BACK
910 /*
911 G4VPhysicalVolume* oldTopVolume =
912 fpTrack->GetTouchableHandle()->GetVolume();
913 fpNavigator->SetNavigatorState(
914 fpITrack->GetTrackingInfo()->GetNavigatorState());
915
916 G4VPhysicalVolume* newTopVolume = fpNavigator->ResetHierarchyAndLocate(
917 fpTrack->GetPosition(), fpTrack->GetMomentumDirection(),
918 *((G4TouchableHistory*) fpTrack->GetTouchableHandle()()));
919
920 // G4VPhysicalVolume* newTopVolume=
921 // fpNavigator->LocateGlobalPointAndSetup(fpTrack->GetPosition(),
922 // &fpTrack->GetMomentumDirection(),
923 // true, false);
924
925 // G4cout << "New Top Volume : " << newTopVolume->GetName() << G4endl;
926
927 if (newTopVolume != oldTopVolume || oldTopVolume->GetRegularStructureId()
928 == 1)
929 {
930 fpState->fTouchableHandle = fpNavigator->CreateTouchableHistory();
931 fpTrack->SetTouchableHandle(fpState->fTouchableHandle);
932 fpTrack->SetNextTouchableHandle(fpState->fTouchableHandle);
933 }
934
935 */
936 //! ADDED BACK
937 //==========================================================================
938 // Only reset navigator state + reset volume hierarchy (internal)
939 // No need to relocate
940 //==========================================================================
941 fpNavigator->SetNavigatorState(fpITrack->GetTrackingInfo()
942 ->GetNavigatorState());
943 }
944 }
945
946 //______________________________________________________________________________
947
948 // ------------------------------------------------------------------------
949 // Compute Interaction Length
950 // ------------------------------------------------------------------------
951 void G4ITStepProcessor::DoDefinePhysicalStepLength()
952 {
953
954 InitDefineStep();
955
956 #ifdef G4VERBOSE
957 // !!!!! Verbose
958 if(fpVerbose != nullptr) fpVerbose->DPSLStarted();
959 #endif
960
961 G4TrackStatus trackStatus = fpTrack->GetTrackStatus();
962
963 if(trackStatus == fStopAndKill)
964 {
965 return;
966 }
967
968 if(trackStatus == fStopButAlive)
969 {
970 fpITrack->GetTrackingInfo()->SetNavigatorState(fpNavigator
971 ->GetNavigatorState());
972 fpNavigator->ResetNavigatorState();
973 return GetAtRestIL();
974 }
975
976 // Find minimum Step length and corresponding time
977 // demanded by active disc./cont. processes
978
979 // ReSet the counter etc.
980 fpState->fPhysicalStep = DBL_MAX; // Initialize by a huge number
981 fPhysIntLength = DBL_MAX; // Initialize by a huge number
982
983 G4double proposedTimeStep = DBL_MAX;
984 G4VProcess* processWithPostStepGivenByTimeStep(nullptr);
985
986 // GPIL for PostStep
987 fPostStepDoItProcTriggered = fpProcessInfo->MAXofPostStepLoops;
988 fPostStepAtTimeDoItProcTriggered = fpProcessInfo->MAXofPostStepLoops;
989
990 // G4cout << "fpProcessInfo->MAXofPostStepLoops : "
991 // << fpProcessInfo->MAXofPostStepLoops
992 // << " mol : " << fpITrack -> GetName()
993 // << " id : " << fpTrack->GetTrackID()
994 // << G4endl;
995
996 for(std::size_t np = 0; np < fpProcessInfo->MAXofPostStepLoops; ++np)
997 {
998 fpCurrentProcess = dynamic_cast<G4VITProcess*>((*fpProcessInfo
999 ->fpPostStepGetPhysIntVector)[(G4int)np]);
1000 if(fpCurrentProcess == nullptr)
1001 {
1002 (fpState->fSelectedPostStepDoItVector)[np] = InActivated;
1003 continue;
1004 } // NULL means the process is inactivated by a user on fly.
1005
1006 fCondition = NotForced;
1007 fpCurrentProcess->SetProcessState(fpTrackingInfo->GetProcessState(fpCurrentProcess
1008 ->GetProcessID()));
1009
1010 // G4cout << "Is going to call : "
1011 // << fpCurrentProcess -> GetProcessName()
1012 // << G4endl;
1013 fPhysIntLength = fpCurrentProcess->PostStepGPIL(*fpTrack,
1014 fpState->fPreviousStepSize,
1015 &fCondition);
1016
1017 #ifdef G4VERBOSE
1018 // !!!!! Verbose
1019 if(fpVerbose != nullptr) fpVerbose->DPSLPostStep();
1020 #endif
1021
1022 fpCurrentProcess->ResetProcessState();
1023 //fpCurrentProcess->SetProcessState(0);
1024
1025 switch(fCondition)
1026 {
1027 case ExclusivelyForced: // Will need special treatment
1028 (fpState->fSelectedPostStepDoItVector)[np] = ExclusivelyForced;
1029 fpState->fStepStatus = fExclusivelyForcedProc;
1030 fpStep->GetPostStepPoint()->SetProcessDefinedStep(fpCurrentProcess);
1031 break;
1032
1033 case Conditionally:
1034 // (fpState->fSelectedPostStepDoItVector)[np] = Conditionally;
1035 G4Exception("G4ITStepProcessor::DefinePhysicalStepLength()",
1036 "ITStepProcessor0008",
1037 FatalException,
1038 "This feature is no more supported");
1039 break;
1040
1041 case Forced:
1042 (fpState->fSelectedPostStepDoItVector)[np] = Forced;
1043 break;
1044
1045 case StronglyForced:
1046 (fpState->fSelectedPostStepDoItVector)[np] = StronglyForced;
1047 break;
1048
1049 default:
1050 (fpState->fSelectedPostStepDoItVector)[np] = InActivated;
1051 break;
1052 }
1053
1054 if(fCondition == ExclusivelyForced)
1055 {
1056 for(std::size_t nrest = np + 1; nrest < fpProcessInfo->MAXofPostStepLoops;
1057 ++nrest)
1058 {
1059 (fpState->fSelectedPostStepDoItVector)[nrest] = InActivated;
1060 }
1061 return; // Please note the 'return' at here !!!
1062 }
1063
1064 if(fPhysIntLength < fpState->fPhysicalStep)
1065 {
1066 // To avoid checking whether the process is actually
1067 // proposing a time step, the returned time steps are
1068 // negative (just for tagging)
1069 if(fpCurrentProcess->ProposesTimeStep())
1070 {
1071 fPhysIntLength *= -1;
1072 if(fPhysIntLength < proposedTimeStep)
1073 {
1074 proposedTimeStep = fPhysIntLength;
1075 fPostStepAtTimeDoItProcTriggered = np;
1076 processWithPostStepGivenByTimeStep = fpCurrentProcess;
1077 }
1078 }
1079 else
1080 {
1081 fpState->fPhysicalStep = fPhysIntLength;
1082 fpState->fStepStatus = fPostStepDoItProc;
1083 fPostStepDoItProcTriggered = G4int(np);
1084 fpStep->GetPostStepPoint()->SetProcessDefinedStep(fpCurrentProcess);
1085 }
1086 }
1087 }
1088
1089 // GPIL for AlongStep
1090 fpState->fProposedSafety = DBL_MAX;
1091 G4double safetyProposedToAndByProcess = fpState->fProposedSafety;
1092
1093 for(std::size_t kp = 0; kp < fpProcessInfo->MAXofAlongStepLoops; ++kp)
1094 {
1095 fpCurrentProcess = dynamic_cast<G4VITProcess*>((*fpProcessInfo
1096 ->fpAlongStepGetPhysIntVector)[(G4int)kp]);
1097 if(fpCurrentProcess == nullptr) continue;
1098 // NULL means the process is inactivated by a user on fly.
1099
1100 fpCurrentProcess->SetProcessState(
1101 fpTrackingInfo->GetProcessState(fpCurrentProcess->GetProcessID()));
1102 fPhysIntLength =
1103 fpCurrentProcess->AlongStepGPIL(*fpTrack,
1104 fpState->fPreviousStepSize,
1105 fpState->fPhysicalStep,
1106 safetyProposedToAndByProcess,
1107 &fGPILSelection);
1108
1109 #ifdef G4VERBOSE
1110 // !!!!! Verbose
1111 if(fpVerbose != nullptr) fpVerbose->DPSLAlongStep();
1112 #endif
1113
1114 if(fPhysIntLength < fpState->fPhysicalStep)
1115 {
1116 fpState->fPhysicalStep = fPhysIntLength;
1117 // Should save PS and TS in IT
1118
1119 // Check if the process wants to be the GPIL winner. For example,
1120 // multi-scattering proposes Step limit, but won't be the winner.
1121 if(fGPILSelection == CandidateForSelection)
1122 {
1123 fpState->fStepStatus = fAlongStepDoItProc;
1124 fpStep->GetPostStepPoint()->SetProcessDefinedStep(fpCurrentProcess);
1125 }
1126
1127 // Transportation is assumed to be the last process in the vector
1128 if(kp == fpProcessInfo->MAXofAlongStepLoops - 1)
1129 {
1130 fpTransportation = dynamic_cast<G4ITTransportation*>(fpCurrentProcess);
1131
1132 if(fpTransportation == nullptr)
1133 {
1134 G4ExceptionDescription exceptionDescription;
1135 exceptionDescription << "No transportation process found ";
1136 G4Exception("G4ITStepProcessor::DoDefinePhysicalStepLength",
1137 "ITStepProcessor0009",
1138 FatalErrorInArgument,
1139 exceptionDescription);
1140 }
1141
1142 fTimeStep = fpTransportation->GetInteractionTimeLeft();
1143
1144 if(fpTrack->GetNextVolume() != nullptr) fpState->fStepStatus = fGeomBoundary;
1145 else fpState->fStepStatus = fWorldBoundary;
1146 }
1147 }
1148 else
1149 {
1150 if(kp == fpProcessInfo->MAXofAlongStepLoops - 1)
1151 {
1152 fpTransportation = dynamic_cast<G4ITTransportation*>(fpCurrentProcess);
1153
1154 if(fpTransportation == nullptr)
1155 {
1156 G4ExceptionDescription exceptionDescription;
1157 exceptionDescription << "No transportation process found ";
1158 G4Exception("G4ITStepProcessor::DoDefinePhysicalStepLength",
1159 "ITStepProcessor0010",
1160 FatalErrorInArgument,
1161 exceptionDescription);
1162 }
1163
1164 fTimeStep = fpTransportation->GetInteractionTimeLeft();
1165 }
1166 }
1167
1168 // Handle PostStep processes sending back time steps rather than space length
1169 if(proposedTimeStep < fTimeStep)
1170 {
1171 if(fPostStepAtTimeDoItProcTriggered < fpProcessInfo->MAXofPostStepLoops)
1172 {
1173 if((fpState->fSelectedPostStepDoItVector)[fPostStepAtTimeDoItProcTriggered] == InActivated)
1174 {
1175 (fpState->fSelectedPostStepDoItVector)[fPostStepAtTimeDoItProcTriggered] =
1176 NotForced;
1177 // (fpState->fSelectedPostStepDoItVector)[fPostStepDoItProcTriggered] = InActivated;
1178
1179 fpState->fStepStatus = fPostStepDoItProc;
1180 fpStep->GetPostStepPoint()->SetProcessDefinedStep(processWithPostStepGivenByTimeStep);
1181
1182 fTimeStep = proposedTimeStep;
1183
1184 fpTransportation->ComputeStep(*fpTrack,
1185 *fpStep,
1186 fTimeStep,
1187 fpState->fPhysicalStep);
1188 }
1189 }
1190 }
1191 else
1192 {
1193 if(fPostStepDoItProcTriggered < fpProcessInfo->MAXofPostStepLoops)
1194 {
1195 if((fpState->fSelectedPostStepDoItVector)[fPostStepDoItProcTriggered] == InActivated)
1196 {
1197 (fpState->fSelectedPostStepDoItVector)[fPostStepDoItProcTriggered] =
1198 NotForced;
1199 }
1200 }
1201 }
1202
1203 // fpCurrentProcess->SetProcessState(0);
1204 fpCurrentProcess->ResetProcessState();
1205
1206 // Make sure to check the safety, even if Step is not limited
1207 // by this process. J. Apostolakis, June 20, 1998
1208 //
1209 if(safetyProposedToAndByProcess < fpState->fProposedSafety)
1210 // proposedSafety keeps the smallest value:
1211 fpState->fProposedSafety = safetyProposedToAndByProcess;
1212 else
1213 // safetyProposedToAndByProcess always proposes a valid safety:
1214 safetyProposedToAndByProcess = fpState->fProposedSafety;
1215
1216 }
1217
1218 fpITrack->GetTrackingInfo()->SetNavigatorState(fpNavigator->GetNavigatorState());
1219 fpNavigator->ResetNavigatorState();
1220 }
1221
1222 //______________________________________________________________________________
1223