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