4 // See the class description in the header file.
6 #include "TG4StepManager.h"
7 #include "TG4GeometryManager.h"
8 #include "TG4PhysicsManager.h"
9 #include "TG4VSensitiveDetector.h"
10 #include "TG4Globals.h"
14 #include <G4SteppingManager.hh>
15 #include <G4UserLimits.hh>
16 #include <G4ParticleTable.hh>
17 #include <G4UImanager.hh>
18 #include <G4AffineTransform.hh>
19 //#include <G4TransportationManager.hh>
20 //#include <G4Navigator.hh>
22 #include <Randomize.hh>
24 #include <TLorentzVector.h>
26 TG4StepManager* TG4StepManager::fgInstance = 0;
28 TG4StepManager::TG4StepManager()
34 TG4Globals::Exception(
35 "TG4StepManager: attempt to create two instances of singleton.");
41 TG4StepManager::TG4StepManager(const TG4StepManager& right) {
43 TG4Globals::Exception(
44 "Attempt to copy TG4StepManager singleton.");
47 TG4StepManager::~TG4StepManager() {
53 TG4StepManager& TG4StepManager::operator=(const TG4StepManager& right)
55 // check assignement to self
56 if (this == &right) return *this;
58 TG4Globals::Exception(
59 "Attempt to assign TG4StepManager singleton.");
66 void TG4StepManager::SetTLorentzVector(G4ThreeVector xyz, G4double t,
67 TLorentzVector& lv) const
69 // Fills TLorentzVector with G4ThreeVector and G4double.
80 void TG4StepManager::StopTrack()
82 // Stops the current track and skips to the next.
83 // ?? do we want to invoke rest processes?
84 // ?? do we want to stop secondaries too?
85 // possible "stop" track status from G4:
86 // fStopButAlive, // Invoke active rest physics processes and
87 // // and kill the current track afterward
88 // fStopAndKill, // Kill the current track
89 // fKillTrackAndSecondaries,
90 // // Kill the current track and also associated
96 fStep->GetTrack()->SetTrackStatus(fStopAndKill);
97 // fStep->GetTrack()->SetTrackStatus(fStopButAlive);
98 // fStep->GetTrack()->SetTrackStatus(fKillTrackAndSecondaries);
101 TG4Globals::Exception("TG4StepManager: Step is not defined.");
105 void TG4StepManager::StopEvent()
107 // Aborts the current event processing.
111 fStep->GetTrack()->SetTrackStatus(fKillTrackAndSecondaries);
112 //StopTrack(); // cannot be used as it keeps secondaries
113 G4UImanager::GetUIpointer()->ApplyCommand("/event/abort");
114 G4UImanager::GetUIpointer()->ApplyCommand("/alStacking/clearStack");
117 TG4Globals::Exception("TG4StepManager: Step is not defined.");
121 void TG4StepManager::Rndm(Float_t* array, const Int_t size) const
123 // Random numbers array of the specified size.
126 G4double* const kpDoubleArray = new G4double[size];
127 RandFlat::shootArray(size,kpDoubleArray);
128 for (G4int i=0; i<size; i++) {
129 array[i] = kpDoubleArray[i];
131 delete [] kpDoubleArray;
134 void TG4StepManager::SetMaxStep(Float_t step)
136 // Maximum step allowed in the current logical volume.
137 // The maximum value is kept for following tracks - is it ok ??
142 G4LogicalVolume* curLogVolume
143 = fStep->GetPreStepPoint()->GetPhysicalVolume()->GetLogicalVolume();
144 G4UserLimits* userLimits
145 = curLogVolume->GetUserLimits();
148 userLimits = new G4UserLimits(step);
149 curLogVolume->SetUserLimits(userLimits);
152 { userLimits->SetMaxAllowedStep(step); }
155 TG4Globals::Exception("TG4StepManager: Step is not defined.");
159 void TG4StepManager::SetMaxNStep(Int_t maxNofSteps)
161 // Not yet implemented.
165 "TG4StepManager::SetMaxNStep(..) is not yet implemented.");
168 void TG4StepManager::SetUserDecay(Int_t pdg)
170 // Not yet implemented.
173 TG4Globals::Exception(
174 "TG4StepManager::SetUserDecay(..) is not yet implemented.");
177 Int_t TG4StepManager::CurrentVolID(Int_t& copyNo) const
179 // Returns the current sensitive detector ID
180 // and the copy number of the current physical volume.
185 G4VPhysicalVolume* physVolume
186 = fStep->GetPreStepPoint()->GetPhysicalVolume();
187 copyNo = physVolume->GetCopyNo();
189 // sensitive detector ID
190 G4VSensitiveDetector* sd
191 = physVolume->GetLogicalVolume()->GetSensitiveDetector();
193 TG4VSensitiveDetector* tsd = dynamic_cast<TG4VSensitiveDetector*>(sd);
197 TG4Globals::Exception(
198 "TG4StepManager::CurrentVol: Unknown sensitive detector type");
203 G4String text = "TG4StepManager::CurrentVol: \n";
204 text = text + " Volume " + physVolume->GetName();
205 text = text + " has not a sensitive detector.";
206 TG4Globals::Exception(text);
211 TG4Globals::Exception("TG4StepManager: Step is not defined.");
216 Int_t TG4StepManager::CurrentVolOffID(Int_t off, Int_t& copyNo) const
218 // Returns the off-th mother's of the current volume
219 // the sensitive detector ID and the copy number.
222 if (off == 0) return CurrentVolID(copyNo);
226 // mother of physical volume may not be set ?!!
228 G4VPhysicalVolume* physVolume
229 = fStep->GetPreStepPoint()->GetPhysicalVolume();
231 G4VPhysicalVolume* mother = 0;
234 mother = physVolume->GetMother();
239 copyNo = mother->GetCopyNo();
241 // sensitive detector ID
242 G4VSensitiveDetector* sd
243 = mother->GetLogicalVolume()->GetSensitiveDetector();
245 TG4VSensitiveDetector* tsd = dynamic_cast<TG4VSensitiveDetector*>(sd);
249 TG4Globals::Exception(
250 "TG4StepManager::CurrentVolOff: Unknown sensitive detector type");
255 G4String text = "TG4StepManager::CurrentVolOff: \n";
256 text = text + " Volume " + mother->GetName();
257 text = text + " has not a sensitive detector.";
258 TG4Globals::Exception(text);
263 G4String text = "TG4StepManager::CurrentVolOff: Volume ";
264 text = text + physVolume->GetName();
265 text = text + " has not defined mother in required level.";
266 TG4Globals::Exception(text);
271 TG4Globals::Exception("TG4StepManager: Step is not defined.");
276 const char* TG4StepManager::CurrentVolName() const
278 // Returns the current physical volume name.
282 G4VPhysicalVolume* physVolume
283 = fStep->GetPreStepPoint()->GetPhysicalVolume();
284 return physVolume->GetName();
287 TG4Globals::Exception("TG4StepManager: Step is not defined.");
292 const char* TG4StepManager::CurrentVolOffName(Int_t off) const
294 // Returns the off-th mother's physical volume name.
298 G4VPhysicalVolume* physVolume
299 = fStep->GetPreStepPoint()->GetPhysicalVolume();
301 G4VPhysicalVolume* mother = 0;
304 mother = physVolume->GetMother();
309 return mother->GetName();
312 TG4Globals::Exception("TG4StepManager::CurrentVolOff: wrong usage.");
317 TG4Globals::Exception("TG4StepManager: Step is not defined.");
322 Int_t TG4StepManager::CurrentMaterial(Float_t &a, Float_t &z, Float_t &dens,
323 Float_t &radl, Float_t &absl) const
325 // Returns the parameters of the current material during transport
326 // the return value is the number of elements in the mixture.
331 G4LogicalVolume* curLogVolume
332 = fStep->GetPreStepPoint()->GetPhysicalVolume()->GetLogicalVolume();
334 = curLogVolume->GetMaterial();
336 // this if may be redundant - check
339 G4int nofElements = material->GetNumberOfElements();
340 TG4GeometryManager* pGeometryManager = TG4GeometryManager::Instance();
341 a = pGeometryManager->GetEffA(material);
342 z = pGeometryManager->GetEffZ(material);
345 dens = material->GetDensity();
346 dens /= TG3Units::MassDensity();
349 radl = material->GetRadlen();
350 radl /= TG3Units::Length();
352 absl = 0.; // this parameter is not defined in Geant4
356 TG4Globals::Exception(
357 "TG4StepManager::CurrentMaterial(..): material is not defined.");
362 TG4Globals::Exception("TG4StepManager: Step is not defined.");
367 void TG4StepManager::Gmtod(Float_t* xm, Float_t* xd, Int_t iflag)
369 // Transforms a position from the world reference frame
370 // to the current volume reference frame.
372 // Geant3 desription:
373 // ==================
374 // Computes coordinates XD (in DRS)
375 // from known coordinates XM in MRS
376 // The local reference system can be initialized by
377 // - the tracking routines and GMTOD used in GUSTEP
378 // - a call to GMEDIA(XM,NUMED)
379 // - a call to GLVOLU(NLEVEL,NAMES,NUMBER,IER)
380 // (inverse routine is GDTOM)
382 // If IFLAG=1 convert coordinates
383 // IFLAG=2 convert direction cosinus
389 G4ThreeVector theGlobalPoint(xm[0],xm[1],xm[2]);
391 //const G4NavigationHistory* history
392 // = fStep->GetPreStepPoint()->GetTouchable()->GetHistory();
393 G4AffineTransform affineTransform
394 = fStep->GetPreStepPoint()->GetTouchable()->GetHistory()
397 G4ThreeVector theLocalPoint;
399 theLocalPoint = affineTransform.TransformPoint(theGlobalPoint);
400 else if ( iflag == 2)
401 theLocalPoint = affineTransform.TransformAxis(theGlobalPoint);
403 TG4Globals::Exception(
404 "TG4StepManager::Gmtod(..,iflag): iflag is not in 1..2");
406 xd[0] = theLocalPoint.x();
407 xd[1] = theLocalPoint.y();
408 xd[2] = theLocalPoint.z();
412 G4ThreeVector direction(0,0,0);
413 G4bool RelativeSearch = true;
414 G4Navigator* theNavigator =
415 G4TransportationManager::GetTransportationManager()->
416 GetNavigatorForTracking();
418 G4VPhysicalVolume* pPhysVol;
420 = LocateGlobalPointAndSetup(theGlobalPoint, &direction, RelativeSearch);
421 //LocateGlobalPointWithinVolume(theGlobalPoint);
424 = theNavigator->GetGlobalToLocalTransform();
428 TG4Globals::Exception("TG4StepManager: Step is not defined.");
432 void TG4StepManager::Gdtom(Float_t* xd, Float_t* xm, Int_t iflag)
434 // Transforms a position from the current volume reference frame
435 // to the world reference frame.
437 // Geant3 desription:
438 // ==================
439 // Computes coordinates XM (Master Reference System
440 // knowing the coordinates XD (Detector Ref System)
441 // The local reference system can be initialized by
442 // - the tracking routines and GDTOM used in GUSTEP
443 // - a call to GSCMED(NLEVEL,NAMES,NUMBER)
444 // (inverse routine is GMTOD)
446 // If IFLAG=1 convert coordinates
447 // IFLAG=2 convert direction cosinus
454 G4ThreeVector theLocalPoint(xd[0],xd[1],xd[2]);
456 G4AffineTransform affineTransform
457 = fStep->GetPreStepPoint()->GetTouchable()->GetHistory()
458 ->GetTopTransform().Inverse();
460 G4ThreeVector theGlobalPoint;
462 theGlobalPoint = affineTransform.TransformPoint(theLocalPoint);
464 theGlobalPoint = affineTransform.TransformAxis(theLocalPoint);
467 "TG4StepManager::Gdtom(...,iflag): iflag is not in 1..2");
469 xm[0] = theGlobalPoint.x();
470 xm[1] = theGlobalPoint.y();
471 xm[2] = theGlobalPoint.z();
474 TG4Globals::Exception("TG4StepManager: Step is not defined.");
478 Float_t TG4StepManager::MaxStep() const
480 // Returns maximum step allowed in the current logical volume
486 G4LogicalVolume* curLogVolume
487 = fStep->GetPreStepPoint()->GetPhysicalVolume()->GetLogicalVolume();
488 G4UserLimits* userLimits
489 = curLogVolume->GetUserLimits();
493 G4String text = "User Limits are not defined for log volume ";
494 text = text + curLogVolume->GetName();
495 TG4Globals::Warning(text);
500 const G4Track& trackRef = *(fStep->GetTrack());
501 maxStep = userLimits->GetMaxAllowedStep(trackRef);
502 maxStep /= TG3Units::Length();
507 TG4Globals::Exception("TG4StepManager: Step is not defined.");
512 Int_t TG4StepManager::GetMaxNStep() const
514 // Not yet implemented.
518 "Method GetMaxNStep is not yet implemented in TG4StepManager.");
522 void TG4StepManager::TrackPosition(TLorentzVector& position) const
524 // Current particle position (in the world reference frame)
525 // and the local time since the current track is created
526 // (position of the PostStepPoint).
531 G4ThreeVector positionVector
532 = fStep->GetPostStepPoint()->GetPosition();
533 positionVector *= 1./(TG3Units::Length());
537 = fStep->GetTrack()->GetLocalTime();
538 time /= TG3Units::Time();
540 SetTLorentzVector(positionVector, time, position);
543 TG4Globals::Exception("TG4StepManager: Step is not defined.");
547 Int_t TG4StepManager::GetMedium() const
549 // Returns the second index of the current material (corresponding to
550 // G3 tracking medium index).
555 G4Material* curMaterial
556 = fStep->GetPreStepPoint()->GetPhysicalVolume()->GetLogicalVolume()->GetMaterial();
559 TG4GeometryManager* pGeometryManager = TG4GeometryManager::Instance();
560 return pGeometryManager->GetMediumId(curMaterial);
563 TG4Globals::Exception("TG4StepManager: Step is not defined.");
568 void TG4StepManager::TrackMomentum(TLorentzVector& momentum) const
570 // Current particle "momentum" (px, py, pz, Etot).
574 G4ThreeVector momentumVector
575 = fStep->GetTrack()->GetMomentum();
578 = fStep->GetTrack()->GetDynamicParticle()->GetTotalEnergy();
579 energy /= TG3Units::Energy();
581 SetTLorentzVector(momentumVector, energy, momentum);
584 TG4Globals::Exception("TG4StepManager: Step is not defined.");
588 void TG4StepManager::TrackVertexPosition(TLorentzVector& position) const
590 // The vertex particle position (in the world reference frame)
591 // and the local time since the current track is created.
596 G4ThreeVector positionVector
597 = fStep->GetTrack()->GetVertexPosition();
598 positionVector *= 1./(TG3Units::Length());
603 = fStep->GetTrack()->GetLocalTime();
604 time /= TG3Units::Time();
606 SetTLorentzVector(positionVector, time, position);
609 TG4Globals::Exception("TG4StepManager: Step is not defined.");
613 void TG4StepManager::TrackVertexMomentum(TLorentzVector& momentum) const
615 // The vertex particle "momentum" (px, py, pz, Ekin)
616 // to do: change Ekin -> Etot
620 G4ThreeVector momentumVector
621 = fStep->GetTrack()->GetVertexMomentumDirection();
624 = fStep->GetTrack()->GetVertexKineticEnergy();
625 energy /= TG3Units::Energy();
627 SetTLorentzVector(momentumVector, energy, momentum);
630 TG4Globals::Exception("TG4StepManager: Step is not defined.");
634 Float_t TG4StepManager::TrackStep() const
636 // Returns the current step length.
641 G4double length = fStep->GetStepLength();
642 length /= TG3Units::Length();
646 TG4Globals::Exception("TG4StepManager: Step is not defined.");
651 Float_t TG4StepManager::TrackLength() const
653 // Returns the length of the current track from its origin.
658 = fStep->GetTrack()->GetTrackLength();
659 length /= TG3Units::Length();
663 TG4Globals::Exception("TG4StepManager: Step is not defined.");
668 Float_t TG4StepManager::TrackTime() const
670 // Returns the local time since the current track is created.
672 // in Geant4: there is also defined proper time as
673 // the proper time of the dynamical particle of the current track.
678 = fStep->GetTrack()->GetLocalTime();
679 time /= TG3Units::Time();
683 TG4Globals::Exception("TG4StepManager: Step is not defined.");
688 Float_t TG4StepManager::Edep() const
690 // Returns total energy deposit in this step.
694 // G4double deltaEnergy = fStep->GetDeltaEnergy();
695 G4double deltaEnergy = fStep->GetTotalEnergyDeposit();
696 deltaEnergy /= TG3Units::Energy();
700 TG4Globals::Exception("TG4StepManager: Step is not defined.");
705 Int_t TG4StepManager::TrackPid() const
707 // Returns the current particle PDG encoding.
711 G4ParticleDefinition* particle
712 = fStep->GetTrack()->GetDynamicParticle()->GetDefinition();
714 // ask TG4PhysicsManager to get PDG encoding
715 // (in order to get PDG from extended TDatabasePDG
716 // in case the standard PDG code is not defined)
717 TG4PhysicsManager* pPhysicsManager = TG4PhysicsManager::Instance();
718 G4int pdgEncoding = pPhysicsManager->GetPDGEncodingFast(particle);
723 TG4Globals::Exception("TG4StepManager: Step is not defined.");
728 Float_t TG4StepManager::TrackCharge() const
730 // Returns the current particle charge.
735 = fStep->GetTrack()->GetDynamicParticle()->GetDefinition()
737 charge /= TG3Units::Charge();
741 TG4Globals::Exception("TG4StepManager: Step is not defined.");
746 Float_t TG4StepManager::TrackMass() const
748 // Returns current particle rest mass.
753 = fStep->GetTrack()->GetDynamicParticle()->GetDefinition()
755 mass /= TG3Units::Mass();
759 TG4Globals::Exception("TG4StepManager: Step is not defined.");
764 Float_t TG4StepManager::Etot() const
766 // Returns total energy of the current particle.
771 = fStep->GetTrack()->GetDynamicParticle()->GetTotalEnergy();
772 energy /= TG3Units::Energy();
776 TG4Globals::Exception("TG4StepManager: Step is not defined.");
781 Bool_t TG4StepManager::IsTrackInside() const
783 // Returns true if particle does not cross geometrical boundary
784 // at both pre-step and post-step points.
787 if ( !(IsTrackEntering()) && !(IsTrackExiting()) )
793 Bool_t TG4StepManager::IsTrackEntering() const
795 // Returns true if particle cross a geometrical boundary
796 // at pre-step point.
801 = fStep->GetPreStepPoint()->GetStepStatus();
802 if (status == fGeomBoundary)
808 TG4Globals::Exception("TG4StepManager: Step is not defined.");
813 Bool_t TG4StepManager::IsTrackExiting() const
815 // Returns true if particle cross a geometrical boundary
816 // at post-step point.
822 = fStep->GetPostStepPoint()->GetStepStatus();
823 if (status == fGeomBoundary)
829 TG4Globals::Exception("TG4StepManager: Step is not defined.");
834 Bool_t TG4StepManager::IsTrackOut() const
836 // Returns true if particle cross the Hall frame
837 // at post-step point.
843 = fStep->GetPostStepPoint()->GetPhysicalVolume()->GetName();
844 if (volName == "pHallFrame")
850 TG4Globals::Exception("TG4StepManager: Step is not defined.");
855 Bool_t TG4StepManager::IsTrackStop() const
857 // Returns true if particle has stopped
858 // or has been killed, suspended or postponed to next event.
860 // Possible track status from G4:
861 // fAlive, // Continue the tracking
862 // fStopButAlive, // Invoke active rest physics processes and
863 // // and kill the current track afterward
864 // fStopAndKill, // Kill the current track
865 // fKillTrackAndSecondaries,
866 // // Kill the current track and also associated
868 // fSuspend, // Suspend the current track
869 // fPostponeToNextEvent
870 // // Postpones the tracking of thecurrent track
871 // // to the next event.
877 = fStep->GetTrack()->GetTrackStatus();
878 if ((status == fStopAndKill) ||
879 (status == fKillTrackAndSecondaries) ||
880 (status == fSuspend) ||
881 (status == fPostponeToNextEvent))
887 TG4Globals::Exception("TG4StepManager: Step is not defined.");
892 Bool_t TG4StepManager::IsTrackDisappeared() const
894 // Returns true if particle has disappeared
895 // (due to any physical process)
896 // or has been killed, suspended or postponed to next event.
902 = fStep->GetTrack()->GetTrackStatus();
903 if ((status == fStopButAlive) ||
904 (status == fKillTrackAndSecondaries) ||
905 (status == fSuspend) ||
906 (status == fPostponeToNextEvent))
912 TG4Globals::Exception("TG4StepManager: Step is not defined.");
917 Bool_t TG4StepManager::IsTrackAlive() const
919 // Returns true if particle continues tracking.
924 = fStep->GetTrack()->GetTrackStatus();
925 if (status == fAlive)
931 TG4Globals::Exception("TG4StepManager: Step is not defined.");
936 Bool_t TG4StepManager::IsNewTrack() const
938 // Returns true when track performs the first step.
941 if (fStep->GetTrack()->GetCurrentStepNumber() == 1)
947 Int_t TG4StepManager::NSecondaries() const
949 // Returns the number of secondary particles generated
950 // in the current step.
953 if (fSteppingManager)
955 G4int nofSecondaries = 0;
956 nofSecondaries += fSteppingManager->GetfN2ndariesAtRestDoIt();
957 nofSecondaries += fSteppingManager->GetfN2ndariesAlongStepDoIt();
958 nofSecondaries += fSteppingManager->GetfN2ndariesPostStepDoIt();
960 return nofSecondaries;
963 TG4Globals::Exception("TG4StepManager: SteppingManager is not defined.");
968 void TG4StepManager::GetSecondary(Int_t index, Int_t& particleId,
969 TLorentzVector& position, TLorentzVector& momentum)
971 // Fills the parameters (particle encoding, position, momentum)
972 // of the generated secondary particle which is specified by index.
973 // !! Check if indexing of secondaries is same !!
976 if (fSteppingManager)
978 G4int nofSecondaries = NSecondaries();
979 G4TrackVector* secondaryTracks = fSteppingManager->GetSecondary();
983 if (index < nofSecondaries)
985 // the index of the first secondary of this step
987 = secondaryTracks->entries() - nofSecondaries;
988 // (the secondaryTracks vector contains secondaries
989 // produced by the track at previous steps, too)
991 = (*secondaryTracks)[startIndex + index];
995 = track->GetDynamicParticle()->GetDefinition()->GetPDGEncoding();
998 G4ThreeVector positionVector = track->GetPosition();
999 positionVector *= 1./(TG3Units::Length());
1000 G4double time = track->GetLocalTime();
1001 time /= TG3Units::Time();
1002 SetTLorentzVector(positionVector, time, position);
1004 // momentum & energy
1005 G4ThreeVector momentumVector = track->GetMomentum();
1006 G4double energy = track->GetDynamicParticle()->GetTotalEnergy();
1007 energy /= TG3Units::Energy();
1008 SetTLorentzVector(momentumVector, energy, momentum);
1011 TG4Globals::Exception(
1012 "TG4StepManager::GetSecondary(): wrong secondary track index.");
1016 TG4Globals::Exception(
1017 "TG4StepManager::GetSecondary(): secondary tracks vector is empty");
1021 TG4Globals::Exception("TG4StepManager: SteppingManager is not defined.");
1025 const char* TG4StepManager::ProdProcess() const
1027 // Returns the name of the process that defined current step
1028 // (and may produce the secondary particles).
1033 const G4VProcess* curProcess
1034 = fStep->GetPostStepPoint()->GetProcessDefinedStep();
1036 G4String g4Name = curProcess->GetProcessName();
1040 TG4Globals::Exception("TG4StepManager: Step is not defined.");