// $Id$
// Category: event
//
+// Author: I.Hrivnacova
+//
+// Class TG4StepManager
+// --------------------
// See the class description in the header file.
#include "TG4StepManager.h"
-#include "TG4GeometryManager.h"
+#include "TG4GeometryServices.h"
+#include "TG4SDServices.h"
+#include "TG4ParticlesManager.h"
#include "TG4PhysicsManager.h"
#include "TG4VSensitiveDetector.h"
#include "TG4Globals.h"
-#include "TG3Units.h"
+#include "TG4G3Units.h"
#include <G4SteppingManager.hh>
#include <G4UserLimits.hh>
#include <G4AffineTransform.hh>
#include <G4TransportationManager.hh>
#include <G4Navigator.hh>
+#include <G4VProcess.hh>
+#include <G4ProcessManager.hh>
+#include <G4ProcessVector.hh>
#include <Randomize.hh>
#include <TLorentzVector.h>
TG4StepManager* TG4StepManager::fgInstance = 0;
+//_____________________________________________________________________________
TG4StepManager::TG4StepManager()
: fTrack(0),
fStep(0),
fgInstance = this;
}
+//_____________________________________________________________________________
TG4StepManager::TG4StepManager(const TG4StepManager& right) {
//
TG4Globals::Exception(
"Attempt to copy TG4StepManager singleton.");
}
+//_____________________________________________________________________________
TG4StepManager::~TG4StepManager() {
//
}
// operators
+//_____________________________________________________________________________
TG4StepManager& TG4StepManager::operator=(const TG4StepManager& right)
{
// check assignement to self
// private methods
+//_____________________________________________________________________________
void TG4StepManager::CheckTrack() const
{
// Gives exception in case the track is not defined.
}
+//_____________________________________________________________________________
void TG4StepManager::CheckStep(const G4String& method) const
{
// Gives exception in case the step is not defined.
}
+//_____________________________________________________________________________
void TG4StepManager::CheckSteppingManager() const
{
// Gives exception in case the step is not defined.
}
+//_____________________________________________________________________________
void TG4StepManager::SetTLorentzVector(G4ThreeVector xyz, G4double t,
TLorentzVector& lv) const
{
lv[3] = t;
}
+//_____________________________________________________________________________
G4VPhysicalVolume* TG4StepManager::GetCurrentOffPhysicalVolume(G4int off) const
{
// Returns the physical volume of the off-th mother's
return mother;
}
-G4int TG4StepManager::GetVolumeID(G4VPhysicalVolume* physVolume) const
-{
-// Returns the sensitive detector ID of the specified
-// physical volume.
-// ---
-
- // sensitive detector ID
- G4VSensitiveDetector* sd
- = physVolume->GetLogicalVolume()->GetSensitiveDetector();
- if (sd) {
- TG4VSensitiveDetector* tsd = dynamic_cast<TG4VSensitiveDetector*>(sd);
- if (tsd)
- return tsd->GetID();
- else {
- TG4Globals::Exception(
- "TG4StepManager::GetVolumeID: Unknown sensitive detector type");
- return 0;
- }
- }
- else {
- G4String text = "TG4StepManager::GetVolumeID: \n";
- text = text + " Volume " + physVolume->GetName();
- text = text + " has not a sensitive detector.";
- //TG4Globals::Exception(text);
- TG4Globals::Warning(text);
- return 0;
- }
-}
-
-
// public methods
+//_____________________________________________________________________________
void TG4StepManager::StopTrack()
{
// Stops the current track and skips to the next.
// fTrack->SetTrackStatus(fKillTrackAndSecondaries);
}
+//_____________________________________________________________________________
void TG4StepManager::StopEvent()
{
// Aborts the current event processing.
G4UImanager::GetUIpointer()->ApplyCommand("/alStacking/clearStack");
}
-void TG4StepManager::Rndm(Float_t* array, const Int_t size) const
-{
-// Random numbers array of the specified size.
-// ---
-
- G4double* const kpDoubleArray = new G4double[size];
- RandFlat::shootArray(size, kpDoubleArray);
- for (G4int i=0; i<size; i++) {
- array[i] = kpDoubleArray[i];
- }
- delete [] kpDoubleArray;
-}
-
+//_____________________________________________________________________________
void TG4StepManager::SetMaxStep(Float_t step)
{
// Maximum step allowed in the current logical volume.
-// The maximum value is kept for following tracks - is it ok ??
// ---
- // check this
- G4LogicalVolume* curLogVolume
- = GetCurrentPhysicalVolume()->GetLogicalVolume();
G4UserLimits* userLimits
- = curLogVolume->GetUserLimits();
- if (userLimits == 0) {
- userLimits = new G4UserLimits(step);
- curLogVolume->SetUserLimits(userLimits);
+ = GetCurrentPhysicalVolume()->GetLogicalVolume()->GetUserLimits();
+
+#ifdef TGEANT4_DEBUG
+ if (!userLimits) {
+ G4String text = "TG4StepManager::SetMaxStep:\n";
+ text = text + " User limits not defined.";
+ TG4Globals::Exception(text);
+ return;
}
- else
- userLimits->SetMaxAllowedStep(step);
+#endif
+
+ // set max step
+ userLimits->SetMaxAllowedStep(step*TG4G3Units::Length());
+
}
+//_____________________________________________________________________________
void TG4StepManager::SetMaxNStep(Int_t maxNofSteps)
{
// Not yet implemented.
"TG4StepManager::SetMaxNStep(..) is not yet implemented.");
}
+//_____________________________________________________________________________
void TG4StepManager::SetUserDecay(Int_t pdg)
{
// Not yet implemented.
"TG4StepManager::SetUserDecay(..) is not yet implemented.");
}
+//_____________________________________________________________________________
G4VPhysicalVolume* TG4StepManager::GetCurrentPhysicalVolume() const
{
// Returns the current physical volume.
return physVolume;
}
+//_____________________________________________________________________________
Int_t TG4StepManager::CurrentVolID(Int_t& copyNo) const
{
// Returns the current sensitive detector ID
copyNo = physVolume->GetCopyNo() + 1;
// sensitive detector ID
- return GetVolumeID(physVolume);
+ TG4SDServices* sdServices = TG4SDServices::Instance();
+ return sdServices->GetVolumeID(physVolume->GetLogicalVolume());
}
+//_____________________________________________________________________________
Int_t TG4StepManager::CurrentVolOffID(Int_t off, Int_t& copyNo) const
{
// Returns the off-th mother's of the current volume
copyNo = mother->GetCopyNo() + 1;
// sensitive detector ID
- return GetVolumeID(mother);
+ TG4SDServices* sdServices = TG4SDServices::Instance();
+ return sdServices->GetVolumeID(mother->GetLogicalVolume());
}
else {
copyNo = 0;
}
}
+//_____________________________________________________________________________
const char* TG4StepManager::CurrentVolName() const
{
// Returns the current physical volume name.
return GetCurrentPhysicalVolume()->GetName();
}
+//_____________________________________________________________________________
const char* TG4StepManager::CurrentVolOffName(Int_t off) const
{
// Returns the off-th mother's physical volume name.
return 0;
}
+//_____________________________________________________________________________
Int_t TG4StepManager::CurrentMaterial(Float_t &a, Float_t &z, Float_t &dens,
Float_t &radl, Float_t &absl) const
{
G4Material* material
= physVolume->GetLogicalVolume()->GetMaterial();
- if (material) {
- G4int nofElements = material->GetNumberOfElements();
- TG4GeometryManager* pGeometryManager = TG4GeometryManager::Instance();
- a = pGeometryManager->GetEffA(material);
- z = pGeometryManager->GetEffZ(material);
-
- // density
- dens = material->GetDensity();
- dens /= TG3Units::MassDensity();
-
- // radiation length
- radl = material->GetRadlen();
- radl /= TG3Units::Length();
-
- absl = 0.; // this parameter is not defined in Geant4
- return nofElements;
- }
- else {
+#ifdef TGEANT4_DEBUG
+ if (!material) {
TG4Globals::Exception(
"TG4StepManager::CurrentMaterial(..): material is not defined.");
return 0;
}
+#endif
+
+ G4int nofElements = material->GetNumberOfElements();
+ TG4GeometryServices* geometryServices = TG4GeometryServices::Instance();
+ a = geometryServices->GetEffA(material);
+ z = geometryServices->GetEffZ(material);
+
+ // density
+ dens = material->GetDensity();
+ dens /= TG4G3Units::MassDensity();
+
+ // radiation length
+ radl = material->GetRadlen();
+ radl /= TG4G3Units::Length();
+
+ absl = 0.; // this parameter is not defined in Geant4
+ return nofElements;
}
+//_____________________________________________________________________________
void TG4StepManager::Gmtod(Float_t* xm, Float_t* xd, Int_t iflag)
{
// Transforms a position from the world reference frame
//
// ---
+#ifdef TGEANT4_DEBUG
+ if (iflag != 1 && iflag != 2) {
+ TG4Globals::Exception(
+ "TG4StepManager::Gmtod(..,iflag): iflag is not in 1..2");
+ return;
+ }
+#endif
+
G4AffineTransform affineTransform;
if (fStepStatus == kVertex) {
G4ThreeVector theGlobalPoint(xm[0],xm[1],xm[2]);
G4ThreeVector theLocalPoint;
- if(iflag == 1)
- theLocalPoint = affineTransform.TransformPoint(theGlobalPoint);
- else if ( iflag == 2)
- theLocalPoint = affineTransform.TransformAxis(theGlobalPoint);
- else
- TG4Globals::Exception(
- "TG4StepManager::Gmtod(..,iflag): iflag is not in 1..2");
+ if (iflag == 1)
+ theLocalPoint = affineTransform.TransformPoint(theGlobalPoint);
+ else
+ // if ( iflag == 2)
+ theLocalPoint = affineTransform.TransformAxis(theGlobalPoint);
xd[0] = theLocalPoint.x();
xd[1] = theLocalPoint.y();
xd[2] = theLocalPoint.z();
}
+//_____________________________________________________________________________
void TG4StepManager::Gdtom(Float_t* xd, Float_t* xm, Int_t iflag)
{
// Transforms a position from the current volume reference frame
xm[2] = theGlobalPoint.z();
}
+//_____________________________________________________________________________
Float_t TG4StepManager::MaxStep() const
{
// Returns maximum step allowed in the current logical volume
else {
const G4Track& trackRef = *(fTrack);
maxStep = userLimits->GetMaxAllowedStep(trackRef);
- maxStep /= TG3Units::Length();
+ maxStep /= TG4G3Units::Length();
return maxStep;
}
}
+//_____________________________________________________________________________
Int_t TG4StepManager::GetMaxNStep() const
{
// Not yet implemented.
return 0;
}
+//_____________________________________________________________________________
void TG4StepManager::TrackPosition(TLorentzVector& position) const
{
// Current particle position (in the world reference frame)
// get position
// check if this is == to PostStepPoint position !!
G4ThreeVector positionVector = fTrack->GetPosition();
- positionVector *= 1./(TG3Units::Length());
+ positionVector *= 1./(TG4G3Units::Length());
// local time
G4double time = fTrack->GetLocalTime();
- time /= TG3Units::Time();
+ time /= TG4G3Units::Time();
SetTLorentzVector(positionVector, time, position);
}
+//_____________________________________________________________________________
Int_t TG4StepManager::GetMedium() const
{
// Returns the second index of the current material (corresponding to
// G3 tracking medium index).
// ---
- // current material
- G4Material* curMaterial
- = GetCurrentPhysicalVolume()->GetLogicalVolume()->GetMaterial();
+ // current logical volume
+ G4LogicalVolume* curLV = GetCurrentPhysicalVolume()->GetLogicalVolume();
// medium index
- TG4GeometryManager* pGeometryManager = TG4GeometryManager::Instance();
- return pGeometryManager->GetMediumId(curMaterial);
+ TG4GeometryServices* geometryServices = TG4GeometryServices::Instance();
+ return geometryServices->GetMediumId(curLV);
}
+//_____________________________________________________________________________
void TG4StepManager::TrackMomentum(TLorentzVector& momentum) const
{
// Current particle "momentum" (px, py, pz, Etot).
#endif
G4ThreeVector momentumVector = fTrack->GetMomentum();
- momentumVector *= 1./(TG3Units::Energy());
+ momentumVector *= 1./(TG4G3Units::Energy());
G4double energy = fTrack->GetDynamicParticle()->GetTotalEnergy();
- energy /= TG3Units::Energy();
+ energy /= TG4G3Units::Energy();
SetTLorentzVector(momentumVector, energy, momentum);
}
+//_____________________________________________________________________________
void TG4StepManager::TrackVertexPosition(TLorentzVector& position) const
{
// The vertex particle position (in the world reference frame)
// position
G4ThreeVector positionVector = fTrack->GetVertexPosition();
- positionVector *= 1./(TG3Units::Length());
+ positionVector *= 1./(TG4G3Units::Length());
// local time
// to be checked
G4double time = fTrack->GetLocalTime();
- time /= TG3Units::Time();
+ time /= TG4G3Units::Time();
SetTLorentzVector(positionVector, time, position);
}
+//_____________________________________________________________________________
void TG4StepManager::TrackVertexMomentum(TLorentzVector& momentum) const
{
// The vertex particle "momentum" (px, py, pz, Ekin)
#endif
G4ThreeVector momentumVector = fTrack->GetVertexMomentumDirection();
- momentumVector *= 1./(TG3Units::Energy());
+ momentumVector *= 1./(TG4G3Units::Energy());
G4double energy = fTrack->GetVertexKineticEnergy();
- energy /= TG3Units::Energy();
+ energy /= TG4G3Units::Energy();
SetTLorentzVector(momentumVector, energy, momentum);
}
+//_____________________________________________________________________________
Float_t TG4StepManager::TrackStep() const
{
// Returns the current step length.
#endif
length = fStep->GetStepLength();
- length /= TG3Units::Length();
+ length /= TG4G3Units::Length();
}
else
length = 0;
return length;
}
+//_____________________________________________________________________________
Float_t TG4StepManager::TrackLength() const
{
// Returns the length of the current track from its origin.
#endif
G4double length = fTrack->GetTrackLength();
- length /= TG3Units::Length();
+ length /= TG4G3Units::Length();
return length;
}
+//_____________________________________________________________________________
Float_t TG4StepManager::TrackTime() const
{
// Returns the local time since the current track is created.
#endif
G4double time = fTrack->GetLocalTime();
- time /= TG3Units::Time();
+ time /= TG4G3Units::Time();
return time;
}
+//_____________________________________________________________________________
Float_t TG4StepManager::Edep() const
{
// Returns total energy deposit in this step.
#endif
energyDeposit = fStep->GetTotalEnergyDeposit();
- energyDeposit /= TG3Units::Energy();
+ energyDeposit /= TG4G3Units::Energy();
}
else
energyDeposit = 0;
return energyDeposit;
}
+//_____________________________________________________________________________
Int_t TG4StepManager::TrackPid() const
{
// Returns the current particle PDG encoding.
G4ParticleDefinition* particle
= fTrack->GetDynamicParticle()->GetDefinition();
- // ask TG4PhysicsManager to get PDG encoding
+ // ask TG4ParticlesManager to get PDG encoding
// (in order to get PDG from extended TDatabasePDG
// in case the standard PDG code is not defined)
- TG4PhysicsManager* pPhysicsManager = TG4PhysicsManager::Instance();
- G4int pdgEncoding = pPhysicsManager->GetPDGEncodingFast(particle);
+ G4int pdgEncoding
+ = TG4ParticlesManager::Instance()->GetPDGEncodingFast(particle);
return pdgEncoding;
}
+//_____________________________________________________________________________
Float_t TG4StepManager::TrackCharge() const
{
// Returns the current particle charge.
G4double charge
= fTrack->GetDynamicParticle()->GetDefinition()
->GetPDGCharge();
- charge /= TG3Units::Charge();
+ charge /= TG4G3Units::Charge();
return charge;
}
+//_____________________________________________________________________________
Float_t TG4StepManager::TrackMass() const
{
// Returns current particle rest mass.
G4double mass
= fTrack->GetDynamicParticle()->GetDefinition()
->GetPDGMass();
- mass /= TG3Units::Mass();
+ mass /= TG4G3Units::Mass();
return mass;
}
+//_____________________________________________________________________________
Float_t TG4StepManager::Etot() const
{
// Returns total energy of the current particle.
G4double energy
= fTrack->GetDynamicParticle()->GetTotalEnergy();
- energy /= TG3Units::Energy();
+ energy /= TG4G3Units::Energy();
return energy;
}
+//_____________________________________________________________________________
Bool_t TG4StepManager::IsTrackInside() const
{
// Returns true if particle does not cross geometrical boundary
return false;
}
+//_____________________________________________________________________________
Bool_t TG4StepManager::IsTrackEntering() const
{
// Returns true if particle cross a geometrical boundary
return false;
}
+//_____________________________________________________________________________
Bool_t TG4StepManager::IsTrackExiting() const
{
// Returns true if particle cross a geometrical boundary.
return false;
}
+//_____________________________________________________________________________
Bool_t TG4StepManager::IsTrackOut() const
{
// Returns true if particle cross the world boundary
return false;
}
+//_____________________________________________________________________________
Bool_t TG4StepManager::IsTrackStop() const
{
// Returns true if particle has stopped
return false;
}
+//_____________________________________________________________________________
Bool_t TG4StepManager::IsTrackDisappeared() const
{
// Returns true if particle has disappeared
return false;
}
+//_____________________________________________________________________________
Bool_t TG4StepManager::IsTrackAlive() const
{
// Returns true if particle continues tracking.
return false;
}
+//_____________________________________________________________________________
Bool_t TG4StepManager::IsNewTrack() const
{
// Returns true when track performs the first step.
return false;
}
+//_____________________________________________________________________________
Int_t TG4StepManager::NSecondaries() const
{
// Returns the number of secondary particles generated
return nofSecondaries;
}
+//_____________________________________________________________________________
void TG4StepManager::GetSecondary(Int_t index, Int_t& particleId,
TLorentzVector& position, TLorentzVector& momentum)
{
G4int nofSecondaries = NSecondaries();
G4TrackVector* secondaryTracks = fSteppingManager->GetSecondary();
- if (secondaryTracks){
- if (index < nofSecondaries) {
-
- // the index of the first secondary of this step
- G4int startIndex
- = secondaryTracks->entries() - nofSecondaries;
- // (the secondaryTracks vector contains secondaries
- // produced by the track at previous steps, too)
- G4Track* track
- = (*secondaryTracks)[startIndex + index];
-
- // particle encoding
- particleId
- = track->GetDynamicParticle()->GetDefinition()->GetPDGEncoding();
-
- // position & time
- G4ThreeVector positionVector = track->GetPosition();
- positionVector *= 1./(TG3Units::Length());
- G4double time = track->GetLocalTime();
- time /= TG3Units::Time();
- SetTLorentzVector(positionVector, time, position);
-
- // momentum & energy
- G4ThreeVector momentumVector = track->GetMomentum();
- G4double energy = track->GetDynamicParticle()->GetTotalEnergy();
- energy /= TG3Units::Energy();
- SetTLorentzVector(momentumVector, energy, momentum);
- }
- else {
- TG4Globals::Exception(
- "TG4StepManager::GetSecondary(): wrong secondary track index.");
- }
- }
- else {
+#ifdef TGEANT4_DEBUG
+ if (!secondaryTracks) {
TG4Globals::Exception(
"TG4StepManager::GetSecondary(): secondary tracks vector is empty");
}
+
+ if (index >= nofSecondaries) {
+ TG4Globals::Exception(
+ "TG4StepManager::GetSecondary(): wrong secondary track index.");
+ }
+#endif
+
+ // the index of the first secondary of this step
+ G4int startIndex
+ = secondaryTracks->size() - nofSecondaries;
+ // (the secondaryTracks vector contains secondaries
+ // produced by the track at previous steps, too)
+ G4Track* track
+ = (*secondaryTracks)[startIndex + index];
+
+ // particle encoding
+ particleId
+ = track->GetDynamicParticle()->GetDefinition()->GetPDGEncoding();
+
+ // position & time
+ G4ThreeVector positionVector = track->GetPosition();
+ positionVector *= 1./(TG4G3Units::Length());
+ G4double time = track->GetLocalTime();
+ time /= TG4G3Units::Time();
+ SetTLorentzVector(positionVector, time, position);
+
+ // momentum & energy
+ G4ThreeVector momentumVector = track->GetMomentum();
+ G4double energy = track->GetDynamicParticle()->GetTotalEnergy();
+ energy /= TG4G3Units::Energy();
+ SetTLorentzVector(momentumVector, energy, momentum);
}
-const char* TG4StepManager::ProdProcess() const
+//_____________________________________________________________________________
+AliMCProcess TG4StepManager::ProdProcess(Int_t isec) const
{
-// Returns the name of the process that defined current step
-// (and may produce the secondary particles).
+// The process that has produced the secondary particles specified
+// with isec index in the current step.
// ---
-
- if (fStepStatus == kVertex) return "NONE";
G4int nofSecondaries = NSecondaries();
- if (nofSecondaries == 0) return "NONE";
-
+ if (fStepStatus == kVertex || !nofSecondaries) return kPNoProcess;
+
#ifdef TGEANT4_DEBUG
CheckStep("ProdProcess");
#endif
- const G4VProcess* curProcess
- = fStep->GetPostStepPoint()->GetProcessDefinedStep();
+ G4TrackVector* secondaryTracks = fSteppingManager->GetSecondary();
+
+#ifdef TGEANT4_DEBUG
+ // should never happen
+ if (!secondaryTracks) {
+ TG4Globals::Exception(
+ "TG4StepManager::ProdProcess(): secondary tracks vector is empty.");
+
+ return kPNoProcess;
+ }
+
+ if (isec >= nofSecondaries) {
+ TG4Globals::Exception(
+ "TG4StepManager::GetSecondary(): wrong secondary track index.");
- G4String g4Name = curProcess->GetProcessName();
- return g4Name;
+ return kPNoProcess;
+ }
+#endif
+
+ // the index of the first secondary of this step
+ G4int startIndex
+ = secondaryTracks->size() - nofSecondaries;
+ // the secondaryTracks vector contains secondaries
+ // produced by the track at previous steps, too
+
+ // the secondary track with specified isec index
+ G4Track* track = (*secondaryTracks)[startIndex + isec];
+
+ const G4VProcess* kpProcess = track->GetCreatorProcess();
+
+ AliMCProcess mcProcess
+ = TG4PhysicsManager::Instance()->GetMCProcess(kpProcess);
+
+ // distinguish kPDeltaRay from kPEnergyLoss
+ if (mcProcess == kPEnergyLoss) mcProcess = kPDeltaRay;
+
+ return mcProcess;
+}
+
+//_____________________________________________________________________________
+Int_t TG4StepManager::StepProcesses(TArrayI& processes) const
+{
+// Fills the array of processes that were active in the current step
+// and returns the number of them.
+// TBD: Distinguish between kPDeltaRay and kPEnergyLoss
+// ---
+
+ if (fStepStatus == kVertex) {
+ G4int nofProcesses = 1;
+ processes.Set(nofProcesses);
+ processes[0] = kPNull;
+ return nofProcesses;
+ }
+
+#ifdef TGEANT4_DEBUG
+ CheckSteppingManager();
+ CheckStep("StepProcesses");
+#endif
+
+ // along step processes
+ G4ProcessVector* processVector
+ = fStep->GetTrack()->GetDefinition()->GetProcessManager()
+ ->GetAlongStepProcessVector();
+ G4int nofAlongStep = processVector->entries();
+
+ // process defined step
+ const G4VProcess* kpLastProcess
+ = fStep->GetPostStepPoint()->GetProcessDefinedStep();
+
+ // set array size
+ processes.Set(nofAlongStep+1);
+ // maximum number of processes:
+ // nofAlongStep (along step) - 1 (transportations) + 1 (post step process)
+ // + possibly 1 (additional process if kPLightScattering)
+ // => nofAlongStep + 1
+
+ // fill array with (nofAlongStep-1) along step processes
+ TG4PhysicsManager* physicsManager = TG4PhysicsManager::Instance();
+ G4int counter = 0;
+ for (G4int i=0; i<nofAlongStep; i++) {
+ G4VProcess* g4Process = (*processVector)[i];
+ // do not fill transportation along step process
+ if (g4Process->GetProcessName() != "Transportation")
+ processes[counter++] = physicsManager->GetMCProcess(g4Process);
+ }
+
+ // fill array with 1 or 2 (if kPLightScattering) last process
+ processes[counter++] = physicsManager->GetMCProcess(kpLastProcess);
+ if (processes[counter-1] == kPLightScattering) {
+ // add reflection/absorption as additional process
+ processes[counter++] = physicsManager->GetOpBoundaryStatus(kpLastProcess);
+ }
+
+ return counter;
}