4 // See the class description in the header file.
6 #include "TG4StepManager.h"
7 #include "TG4GeometryServices.h"
8 #include "TG4PhysicsManager.h"
9 #include "TG4VSensitiveDetector.h"
10 #include "TG4Limits.h"
11 #include "TG4Globals.h"
14 #include <G4SteppingManager.hh>
15 #include <G4UserLimits.hh>
16 #include <G4UImanager.hh>
17 #include <G4AffineTransform.hh>
18 #include <G4TransportationManager.hh>
19 #include <G4Navigator.hh>
20 #include <G4VProcess.hh>
21 #include <G4ProcessManager.hh>
22 #include <G4ProcessVector.hh>
24 #include <Randomize.hh>
25 #include <TLorentzVector.h>
27 TG4StepManager* TG4StepManager::fgInstance = 0;
29 TG4StepManager::TG4StepManager()
32 fStepStatus(kNormalStep),
37 TG4Globals::Exception(
38 "TG4StepManager: attempt to create two instances of singleton.");
44 TG4StepManager::TG4StepManager(const TG4StepManager& right) {
46 TG4Globals::Exception(
47 "Attempt to copy TG4StepManager singleton.");
50 TG4StepManager::~TG4StepManager() {
56 TG4StepManager& TG4StepManager::operator=(const TG4StepManager& right)
58 // check assignement to self
59 if (this == &right) return *this;
61 TG4Globals::Exception(
62 "Attempt to assign TG4StepManager singleton.");
69 void TG4StepManager::CheckTrack() const
71 // Gives exception in case the track is not defined.
75 TG4Globals::Exception("TG4StepManager: Track is not defined.");
79 void TG4StepManager::CheckStep(const G4String& method) const
81 // Gives exception in case the step is not defined.
85 G4String text = "TG4StepManager::";
86 text = text + method + ": Step is not defined.";
87 TG4Globals::Exception(text);
92 void TG4StepManager::CheckSteppingManager() const
94 // Gives exception in case the step is not defined.
97 if (!fSteppingManager)
98 TG4Globals::Exception("TG4StepManager: Stepping manager is not defined.");
102 void TG4StepManager::SetTLorentzVector(G4ThreeVector xyz, G4double t,
103 TLorentzVector& lv) const
105 // Fills TLorentzVector with G4ThreeVector and G4double.
114 G4VPhysicalVolume* TG4StepManager::GetCurrentOffPhysicalVolume(G4int off) const
116 // Returns the physical volume of the off-th mother's
117 // of the current volume.
120 G4VPhysicalVolume* physVolume = GetCurrentPhysicalVolume();
122 G4VPhysicalVolume* mother = physVolume;
126 if (mother) mother = mother->GetMother();
131 G4String text = "TG4StepManager::CurrentVolOff: \n";
132 text = text + " Volume ";
133 text = text + physVolume->GetName();
134 text = text + " has not defined mother in the required level.";
135 TG4Globals::Warning(text);
143 void TG4StepManager::StopTrack()
145 // Stops the current track and skips to the next.
146 // ?? do we want to invoke rest processes?
147 // ?? do we want to stop secondaries too?
148 // possible "stop" track status from G4:
149 // fStopButAlive, // Invoke active rest physics processes and
150 // // and kill the current track afterward
151 // fStopAndKill, // Kill the current track
152 // fKillTrackAndSecondaries,
153 // // Kill the current track and also associated
161 fTrack->SetTrackStatus(fStopAndKill);
162 // fTrack->SetTrackStatus(fStopButAlive);
163 // fTrack->SetTrackStatus(fKillTrackAndSecondaries);
166 void TG4StepManager::StopEvent()
168 // Aborts the current event processing.
175 fTrack->SetTrackStatus(fKillTrackAndSecondaries);
176 //StopTrack(); // cannot be used as it keeps secondaries
177 G4UImanager::GetUIpointer()->ApplyCommand("/event/abort");
178 G4UImanager::GetUIpointer()->ApplyCommand("/alStacking/clearStack");
181 void TG4StepManager::SetMaxStep(Float_t step)
183 // Maximum step allowed in the current logical volume.
186 G4LogicalVolume* curLogVolume
187 = GetCurrentPhysicalVolume()->GetLogicalVolume();
188 G4UserLimits* userLimits
189 = curLogVolume->GetUserLimits();
191 if (userLimits == 0) {
193 userLimits = new TG4Limits();
195 // set limits to all logical volumes
196 // corresponding to the current "G3" volume
197 TG4GeometryServices* geometryServices = TG4GeometryServices::Instance();
198 G4int nofLV = geometryServices->SetUserLimits(userLimits, curLogVolume);
202 userLimits->SetMaxAllowedStep(step*TG3Units::Length());
205 void TG4StepManager::SetMaxNStep(Int_t maxNofSteps)
207 // Not yet implemented.
211 "TG4StepManager::SetMaxNStep(..) is not yet implemented.");
214 void TG4StepManager::SetUserDecay(Int_t pdg)
216 // Not yet implemented.
219 TG4Globals::Exception(
220 "TG4StepManager::SetUserDecay(..) is not yet implemented.");
223 G4VPhysicalVolume* TG4StepManager::GetCurrentPhysicalVolume() const
225 // Returns the current physical volume.
226 // According to fStepStatus the volume from track vertex,
227 // pre step point or post step point is returned.
230 G4VPhysicalVolume* physVolume;
231 if (fStepStatus == kNormalStep) {
234 CheckStep("GetCurrentPhysicalVolume");
237 physVolume = fStep->GetPreStepPoint()->GetPhysicalVolume();
239 else if (fStepStatus == kBoundary) {
242 CheckStep("GetCurrentPhysicalVolume");
245 physVolume = fStep->GetPostStepPoint()->GetPhysicalVolume();
253 G4ThreeVector position = fTrack->GetPosition();
254 G4Navigator* navigator =
255 G4TransportationManager::GetTransportationManager()->
256 GetNavigatorForTracking();
258 = navigator->LocateGlobalPointAndSetup(position);
264 Int_t TG4StepManager::CurrentVolID(Int_t& copyNo) const
266 // Returns the current sensitive detector ID
267 // and the copy number of the current physical volume.
270 G4VPhysicalVolume* physVolume = GetCurrentPhysicalVolume();
271 copyNo = physVolume->GetCopyNo() + 1;
273 // sensitive detector ID
274 TG4GeometryServices* geometryServices = TG4GeometryServices::Instance();
275 return geometryServices->GetVolumeID(physVolume->GetLogicalVolume());
278 Int_t TG4StepManager::CurrentVolOffID(Int_t off, Int_t& copyNo) const
280 // Returns the off-th mother's of the current volume
281 // the sensitive detector ID and the copy number.
284 if (off == 0) return CurrentVolID(copyNo);
286 G4VPhysicalVolume* mother = GetCurrentOffPhysicalVolume(off);
289 copyNo = mother->GetCopyNo() + 1;
291 // sensitive detector ID
292 TG4GeometryServices* geometryServices = TG4GeometryServices::Instance();
293 return geometryServices->GetVolumeID(mother->GetLogicalVolume());
301 const char* TG4StepManager::CurrentVolName() const
303 // Returns the current physical volume name.
306 return GetCurrentPhysicalVolume()->GetName();
309 const char* TG4StepManager::CurrentVolOffName(Int_t off) const
311 // Returns the off-th mother's physical volume name.
314 if (off == 0) return CurrentVolName();
316 G4VPhysicalVolume* mother = GetCurrentOffPhysicalVolume(off);
319 return mother->GetName();
324 Int_t TG4StepManager::CurrentMaterial(Float_t &a, Float_t &z, Float_t &dens,
325 Float_t &radl, Float_t &absl) const
327 // Returns the parameters of the current material during transport
328 // the return value is the number of elements in the mixture.
331 G4VPhysicalVolume* physVolume = GetCurrentPhysicalVolume();
334 = physVolume->GetLogicalVolume()->GetMaterial();
337 G4int nofElements = material->GetNumberOfElements();
338 TG4GeometryServices* geometryServices = TG4GeometryServices::Instance();
339 a = geometryServices->GetEffA(material);
340 z = geometryServices->GetEffZ(material);
343 dens = material->GetDensity();
344 dens /= TG3Units::MassDensity();
347 radl = material->GetRadlen();
348 radl /= TG3Units::Length();
350 absl = 0.; // this parameter is not defined in Geant4
354 TG4Globals::Exception(
355 "TG4StepManager::CurrentMaterial(..): material is not defined.");
360 void TG4StepManager::Gmtod(Float_t* xm, Float_t* xd, Int_t iflag)
362 // Transforms a position from the world reference frame
363 // to the current volume reference frame.
365 // Geant3 desription:
366 // ==================
367 // Computes coordinates XD (in DRS)
368 // from known coordinates XM in MRS
369 // The local reference system can be initialized by
370 // - the tracking routines and GMTOD used in GUSTEP
371 // - a call to GMEDIA(XM,NUMED)
372 // - a call to GLVOLU(NLEVEL,NAMES,NUMBER,IER)
373 // (inverse routine is GDTOM)
375 // If IFLAG=1 convert coordinates
376 // IFLAG=2 convert direction cosinus
380 G4AffineTransform affineTransform;
382 if (fStepStatus == kVertex) {
383 G4Navigator* navigator =
384 G4TransportationManager::GetTransportationManager()->
385 GetNavigatorForTracking();
387 affineTransform = navigator->GetGlobalToLocalTransform();
396 = fStep->GetPreStepPoint()->GetTouchable()->GetHistory()
400 G4ThreeVector theGlobalPoint(xm[0],xm[1],xm[2]);
401 G4ThreeVector theLocalPoint;
403 theLocalPoint = affineTransform.TransformPoint(theGlobalPoint);
404 else if ( iflag == 2)
405 theLocalPoint = affineTransform.TransformAxis(theGlobalPoint);
407 TG4Globals::Exception(
408 "TG4StepManager::Gmtod(..,iflag): iflag is not in 1..2");
410 xd[0] = theLocalPoint.x();
411 xd[1] = theLocalPoint.y();
412 xd[2] = theLocalPoint.z();
415 void TG4StepManager::Gdtom(Float_t* xd, Float_t* xm, Int_t iflag)
417 // Transforms a position from the current volume reference frame
418 // to the world reference frame.
420 // Geant3 desription:
421 // ==================
422 // Computes coordinates XM (Master Reference System
423 // knowing the coordinates XD (Detector Ref System)
424 // The local reference system can be initialized by
425 // - the tracking routines and GDTOM used in GUSTEP
426 // - a call to GSCMED(NLEVEL,NAMES,NUMBER)
427 // (inverse routine is GMTOD)
429 // If IFLAG=1 convert coordinates
430 // IFLAG=2 convert direction cosinus
434 G4AffineTransform affineTransform;
436 if (fStepStatus == kVertex) {
437 G4Navigator* navigator =
438 G4TransportationManager::GetTransportationManager()->
439 GetNavigatorForTracking();
441 affineTransform = navigator->GetLocalToGlobalTransform();
452 = fStep->GetPreStepPoint()->GetTouchable()->GetHistory()
453 ->GetTopTransform().Inverse();
456 G4ThreeVector theLocalPoint(xd[0],xd[1],xd[2]);
457 G4ThreeVector theGlobalPoint;
459 theGlobalPoint = affineTransform.TransformPoint(theLocalPoint);
461 theGlobalPoint = affineTransform.TransformAxis(theLocalPoint);
464 "TG4StepManager::Gdtom(...,iflag): iflag is not in 1..2");
466 xm[0] = theGlobalPoint.x();
467 xm[1] = theGlobalPoint.y();
468 xm[2] = theGlobalPoint.z();
471 Float_t TG4StepManager::MaxStep() const
473 // Returns maximum step allowed in the current logical volume
477 G4LogicalVolume* curLogVolume
478 = GetCurrentPhysicalVolume()->GetLogicalVolume();
481 G4UserLimits* userLimits
482 = curLogVolume->GetUserLimits();
485 if (userLimits == 0) {
486 G4String text = "User Limits are not defined for log volume ";
487 text = text + curLogVolume->GetName();
488 TG4Globals::Warning(text);
492 const G4Track& trackRef = *(fTrack);
493 maxStep = userLimits->GetMaxAllowedStep(trackRef);
494 maxStep /= TG3Units::Length();
499 Int_t TG4StepManager::GetMaxNStep() const
501 // Not yet implemented.
505 "Method GetMaxNStep is not yet implemented in TG4StepManager.");
509 void TG4StepManager::TrackPosition(TLorentzVector& position) const
511 // Current particle position (in the world reference frame)
512 // and the local time since the current track is created
513 // (position of the PostStepPoint).
521 // check if this is == to PostStepPoint position !!
522 G4ThreeVector positionVector = fTrack->GetPosition();
523 positionVector *= 1./(TG3Units::Length());
526 G4double time = fTrack->GetLocalTime();
527 time /= TG3Units::Time();
529 SetTLorentzVector(positionVector, time, position);
532 Int_t TG4StepManager::GetMedium() const
534 // Returns the second index of the current material (corresponding to
535 // G3 tracking medium index).
539 G4Material* curMaterial
540 = GetCurrentPhysicalVolume()->GetLogicalVolume()->GetMaterial();
543 TG4GeometryServices* geometryServices = TG4GeometryServices::Instance();
544 return geometryServices->GetMediumId(curMaterial);
547 void TG4StepManager::TrackMomentum(TLorentzVector& momentum) const
549 // Current particle "momentum" (px, py, pz, Etot).
556 G4ThreeVector momentumVector = fTrack->GetMomentum();
557 momentumVector *= 1./(TG3Units::Energy());
559 G4double energy = fTrack->GetDynamicParticle()->GetTotalEnergy();
560 energy /= TG3Units::Energy();
562 SetTLorentzVector(momentumVector, energy, momentum);
565 void TG4StepManager::TrackVertexPosition(TLorentzVector& position) const
567 // The vertex particle position (in the world reference frame)
568 // and the local time since the current track is created.
576 G4ThreeVector positionVector = fTrack->GetVertexPosition();
577 positionVector *= 1./(TG3Units::Length());
581 G4double time = fTrack->GetLocalTime();
582 time /= TG3Units::Time();
584 SetTLorentzVector(positionVector, time, position);
587 void TG4StepManager::TrackVertexMomentum(TLorentzVector& momentum) const
589 // The vertex particle "momentum" (px, py, pz, Ekin)
590 // to do: change Ekin -> Etot
597 G4ThreeVector momentumVector = fTrack->GetVertexMomentumDirection();
598 momentumVector *= 1./(TG3Units::Energy());
600 G4double energy = fTrack->GetVertexKineticEnergy();
601 energy /= TG3Units::Energy();
603 SetTLorentzVector(momentumVector, energy, momentum);
606 Float_t TG4StepManager::TrackStep() const
608 // Returns the current step length.
612 if (fStepStatus == kNormalStep) {
615 CheckStep("TrackStep");
618 length = fStep->GetStepLength();
619 length /= TG3Units::Length();
627 Float_t TG4StepManager::TrackLength() const
629 // Returns the length of the current track from its origin.
636 G4double length = fTrack->GetTrackLength();
637 length /= TG3Units::Length();
641 Float_t TG4StepManager::TrackTime() const
643 // Returns the local time since the current track is created.
645 // in Geant4: there is also defined proper time as
646 // the proper time of the dynamical particle of the current track.
653 G4double time = fTrack->GetLocalTime();
654 time /= TG3Units::Time();
658 Float_t TG4StepManager::Edep() const
660 // Returns total energy deposit in this step.
663 G4double energyDeposit;
664 if (fStepStatus == kNormalStep) {
670 energyDeposit = fStep->GetTotalEnergyDeposit();
671 energyDeposit /= TG3Units::Energy();
676 return energyDeposit;
679 Int_t TG4StepManager::TrackPid() const
681 // Returns the current particle PDG encoding.
688 G4ParticleDefinition* particle
689 = fTrack->GetDynamicParticle()->GetDefinition();
691 // ask TG4PhysicsManager to get PDG encoding
692 // (in order to get PDG from extended TDatabasePDG
693 // in case the standard PDG code is not defined)
694 TG4PhysicsManager* pPhysicsManager = TG4PhysicsManager::Instance();
695 G4int pdgEncoding = pPhysicsManager->GetPDGEncodingFast(particle);
700 Float_t TG4StepManager::TrackCharge() const
702 // Returns the current particle charge.
710 = fTrack->GetDynamicParticle()->GetDefinition()
712 charge /= TG3Units::Charge();
716 Float_t TG4StepManager::TrackMass() const
718 // Returns current particle rest mass.
726 = fTrack->GetDynamicParticle()->GetDefinition()
728 mass /= TG3Units::Mass();
732 Float_t TG4StepManager::Etot() const
734 // Returns total energy of the current particle.
742 = fTrack->GetDynamicParticle()->GetTotalEnergy();
743 energy /= TG3Units::Energy();
747 Bool_t TG4StepManager::IsTrackInside() const
749 // Returns true if particle does not cross geometrical boundary
750 // and is not in vertex.
753 if (fStepStatus == kNormalStep && !(IsTrackExiting()) ) {
754 // track is always inside during a normal step
761 Bool_t TG4StepManager::IsTrackEntering() const
763 // Returns true if particle cross a geometrical boundary
767 if (fStepStatus != kNormalStep) {
768 // track is entering during a vertex or boundary step
775 Bool_t TG4StepManager::IsTrackExiting() const
777 // Returns true if particle cross a geometrical boundary.
780 if (fStepStatus == kNormalStep) {
783 CheckStep("IsTrackExiting");
786 if (fStep->GetPostStepPoint()->GetStepStatus() == fGeomBoundary)
793 Bool_t TG4StepManager::IsTrackOut() const
795 // Returns true if particle cross the world boundary
796 // at post-step point.
799 if (fStepStatus == kVertex) return false;
802 CheckStep("IsTrackCut");
807 = fStep->GetPostStepPoint()->GetStepStatus();
808 if (status == fWorldBoundary)
814 Bool_t TG4StepManager::IsTrackStop() const
816 // Returns true if particle has stopped
817 // or has been killed, suspended or postponed to next event.
819 // Possible track status from G4:
820 // fAlive, // Continue the tracking
821 // fStopButAlive, // Invoke active rest physics processes and
822 // // and kill the current track afterward
823 // fStopAndKill, // Kill the current track
824 // fKillTrackAndSecondaries,
825 // // Kill the current track and also associated
827 // fSuspend, // Suspend the current track
828 // fPostponeToNextEvent
829 // // Postpones the tracking of thecurrent track
830 // // to the next event.
839 = fTrack->GetTrackStatus();
840 if ((status == fStopAndKill) ||
841 (status == fKillTrackAndSecondaries) ||
842 (status == fSuspend) ||
843 (status == fPostponeToNextEvent)) {
850 Bool_t TG4StepManager::IsTrackDisappeared() const
852 // Returns true if particle has disappeared
853 // (due to any physical process)
854 // or has been killed, suspended or postponed to next event.
863 = fTrack->GetTrackStatus();
864 if ((status == fStopButAlive) ||
865 (status == fKillTrackAndSecondaries) ||
866 (status == fSuspend) ||
867 (status == fPostponeToNextEvent)) {
874 Bool_t TG4StepManager::IsTrackAlive() const
876 // Returns true if particle continues tracking.
884 = fTrack->GetTrackStatus();
885 if (status == fAlive)
891 Bool_t TG4StepManager::IsNewTrack() const
893 // Returns true when track performs the first step.
896 if (fStepStatus == kVertex)
902 Int_t TG4StepManager::NSecondaries() const
904 // Returns the number of secondary particles generated
905 // in the current step.
909 CheckSteppingManager();
912 G4int nofSecondaries = 0;
913 nofSecondaries += fSteppingManager->GetfN2ndariesAtRestDoIt();
914 nofSecondaries += fSteppingManager->GetfN2ndariesAlongStepDoIt();
915 nofSecondaries += fSteppingManager->GetfN2ndariesPostStepDoIt();
917 return nofSecondaries;
920 void TG4StepManager::GetSecondary(Int_t index, Int_t& particleId,
921 TLorentzVector& position, TLorentzVector& momentum)
923 // Fills the parameters (particle encoding, position, momentum)
924 // of the generated secondary particle which is specified by index.
925 // !! Check if indexing of secondaries is same !!
929 CheckSteppingManager();
932 G4int nofSecondaries = NSecondaries();
933 G4TrackVector* secondaryTracks = fSteppingManager->GetSecondary();
935 if (secondaryTracks){
936 if (index < nofSecondaries) {
938 // the index of the first secondary of this step
940 = secondaryTracks->entries() - nofSecondaries;
941 // (the secondaryTracks vector contains secondaries
942 // produced by the track at previous steps, too)
944 = (*secondaryTracks)[startIndex + index];
948 = track->GetDynamicParticle()->GetDefinition()->GetPDGEncoding();
951 G4ThreeVector positionVector = track->GetPosition();
952 positionVector *= 1./(TG3Units::Length());
953 G4double time = track->GetLocalTime();
954 time /= TG3Units::Time();
955 SetTLorentzVector(positionVector, time, position);
958 G4ThreeVector momentumVector = track->GetMomentum();
959 G4double energy = track->GetDynamicParticle()->GetTotalEnergy();
960 energy /= TG3Units::Energy();
961 SetTLorentzVector(momentumVector, energy, momentum);
964 TG4Globals::Exception(
965 "TG4StepManager::GetSecondary(): wrong secondary track index.");
969 TG4Globals::Exception(
970 "TG4StepManager::GetSecondary(): secondary tracks vector is empty");
974 AliMCProcess TG4StepManager::ProdProcess(Int_t isec) const
976 // The process that has produced the secondary particles specified
977 // with isec index in the current step.
980 G4int nofSecondaries = NSecondaries();
981 if (fStepStatus == kVertex || !nofSecondaries) return kPNoProcess;
984 CheckStep("ProdProcess");
987 G4TrackVector* secondaryTracks = fSteppingManager->GetSecondary();
989 // should never happen
990 if (!secondaryTracks) {
991 TG4Globals::Exception(
992 "TG4StepManager::ProdProcess(): secondary tracks vector is empty.");
997 if (isec < nofSecondaries) {
999 // the index of the first secondary of this step
1001 = secondaryTracks->entries() - nofSecondaries;
1002 // the secondaryTracks vector contains secondaries
1003 // produced by the track at previous steps, too
1005 // the secondary track with specified isec index
1006 G4Track* track = (*secondaryTracks)[startIndex + isec];
1008 const G4VProcess* kpProcess = track->GetCreatorProcess();
1010 TG4PhysicsManager* pPhysicsManager = TG4PhysicsManager::Instance();
1011 AliMCProcess mcProcess = pPhysicsManager->GetMCProcess(kpProcess);
1013 // distinguish kPDeltaRay from kPEnergyLoss
1014 if (mcProcess == kPEnergyLoss) mcProcess = kPDeltaRay;
1019 TG4Globals::Exception(
1020 "TG4StepManager::GetSecondary(): wrong secondary track index.");
1027 Int_t TG4StepManager::StepProcesses(TArrayI &proc) const
1029 // Fills the array of processes that were active in the current step
1030 // and returns the number of them.
1031 // TBD: Distinguish between kPDeltaRay and kPEnergyLoss
1034 if (fStepStatus == kVertex) {
1035 G4cout << "kVertex" << G4endl;
1036 G4int nofProcesses = 1;
1037 proc.Set(nofProcesses);
1039 return nofProcesses;
1042 #ifdef TGEANT4_DEBUG
1043 CheckSteppingManager();
1044 CheckStep("StepProcesses");
1047 // along step processes
1048 G4ProcessManager* processManager
1049 = fStep->GetTrack()->GetDefinition()->GetProcessManager();
1050 G4ProcessVector* alongStepProcessVector
1051 = processManager->GetAlongStepProcessVector();
1052 G4int nofProcesses = alongStepProcessVector->entries();
1054 // process defined step
1055 const G4VProcess* kpLastProcess
1056 = fStep->GetPostStepPoint()->GetProcessDefinedStep();
1058 // fill the array of processes
1059 proc.Set(nofProcesses);
1060 TG4PhysicsManager* physicsManager = TG4PhysicsManager::Instance();
1062 for (i=0; i<nofProcesses-1; i++) {
1063 G4VProcess* g4Process = (*alongStepProcessVector)[i];
1064 // do not fill transportation along step process
1065 if (g4Process->GetProcessName() != "Transportation") {
1066 physicsManager->GetMCProcess(g4Process);
1067 proc[i] = physicsManager->GetMCProcess(g4Process);
1070 proc[nofProcesses-1] = physicsManager->GetMCProcess(kpLastProcess);
1072 return nofProcesses;