// See the class description in the header file.
#include "TG4StepManager.h"
-#include "TG4GeometryManager.h"
+#include "TG4GeometryServices.h"
+#include "TG4ParticlesManager.h"
#include "TG4PhysicsManager.h"
#include "TG4VSensitiveDetector.h"
+#include "TG4Limits.h"
#include "TG4Globals.h"
-#include "TG3Units.h"
+#include "TG4G3Units.h"
#include <G4SteppingManager.hh>
#include <G4UserLimits.hh>
-#include <G4ParticleTable.hh>
#include <G4UImanager.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;
}
-void TG4StepManager::CheckStep() const
+void TG4StepManager::CheckStep(const G4String& method) const
{
// Gives exception in case the step is not defined.
// ---
- if (!fStep)
- TG4Globals::Exception("TG4StepManager: Step is not defined.");
+ if (!fStep) {
+ G4String text = "TG4StepManager::";
+ text = text + method + ": Step is not defined.";
+ TG4Globals::Exception(text);
+ }
}
G4VPhysicalVolume* physVolume = GetCurrentPhysicalVolume();
- G4VPhysicalVolume* mother = physVolume;
+ G4VPhysicalVolume* mother = physVolume;
+
Int_t level = off;
while (level > 0) {
if (mother) mother = mother->GetMother();
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::CurrentVol: \n";
- text = text + " Volume " + physVolume->GetName();
- text = text + " has not a sensitive detector.";
- TG4Globals::Exception(text);
- return 0;
- }
-}
-
-
// public methods
void TG4StepManager::StopTrack()
// // secondaries.
// ---
+#ifdef TGEANT4_DEBUG
CheckTrack();
+#endif
fTrack->SetTrackStatus(fStopAndKill);
// fTrack->SetTrackStatus(fStopButAlive);
// Aborts the current event processing.
// ---
+#ifdef TGEANT4_DEBUG
CheckTrack();
+#endif
fTrack->SetTrackStatus(fKillTrackAndSecondaries);
//StopTrack(); // cannot be used as it keeps secondaries
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);
+ // create new limits
+ userLimits = new TG4Limits();
+
+ // set limits to all logical volumes
+ // corresponding to the current "G3" volume
+ TG4GeometryServices* geometryServices = TG4GeometryServices::Instance();
+ G4int nofLV = geometryServices->SetUserLimits(userLimits, curLogVolume);
}
- else
- userLimits->SetMaxAllowedStep(step);
+
+ // set max step
+ userLimits->SetMaxAllowedStep(step*TG4G3Units::Length());
}
void TG4StepManager::SetMaxNStep(Int_t maxNofSteps)
G4VPhysicalVolume* physVolume;
if (fStepStatus == kNormalStep) {
- CheckStep();
+
+#ifdef TGEANT4_DEBUG
+ CheckStep("GetCurrentPhysicalVolume");
+#endif
+
physVolume = fStep->GetPreStepPoint()->GetPhysicalVolume();
}
else if (fStepStatus == kBoundary) {
- CheckStep();
+
+#ifdef TGEANT4_DEBUG
+ CheckStep("GetCurrentPhysicalVolume");
+#endif
+
physVolume = fStep->GetPostStepPoint()->GetPhysicalVolume();
}
else {
+
+#ifdef TGEANT4_DEBUG
CheckTrack();
+#endif
+
G4ThreeVector position = fTrack->GetPosition();
G4Navigator* navigator =
G4TransportationManager::GetTransportationManager()->
// ---
G4VPhysicalVolume* physVolume = GetCurrentPhysicalVolume();
- copyNo = physVolume->GetCopyNo();
+ copyNo = physVolume->GetCopyNo() + 1;
// sensitive detector ID
- return GetVolumeID(physVolume);
+ TG4GeometryServices* geometryServices = TG4GeometryServices::Instance();
+ return geometryServices->GetVolumeID(physVolume->GetLogicalVolume());
}
Int_t TG4StepManager::CurrentVolOffID(Int_t off, Int_t& copyNo) const
G4VPhysicalVolume* mother = GetCurrentOffPhysicalVolume(off);
if (mother) {
- copyNo = mother->GetCopyNo();
+ copyNo = mother->GetCopyNo() + 1;
// sensitive detector ID
- return GetVolumeID(mother);
+ TG4GeometryServices* geometryServices = TG4GeometryServices::Instance();
+ return geometryServices->GetVolumeID(mother->GetLogicalVolume());
}
else {
copyNo = 0;
G4Material* material
= physVolume->GetLogicalVolume()->GetMaterial();
- // this if may be redundant - check
- if (material)
- {
+ if (material) {
G4int nofElements = material->GetNumberOfElements();
- TG4GeometryManager* pGeometryManager = TG4GeometryManager::Instance();
- a = pGeometryManager->GetEffA(material);
- z = pGeometryManager->GetEffZ(material);
+ TG4GeometryServices* geometryServices = TG4GeometryServices::Instance();
+ a = geometryServices->GetEffA(material);
+ z = geometryServices->GetEffZ(material);
// density
dens = material->GetDensity();
- dens /= TG3Units::MassDensity();
+ dens /= TG4G3Units::MassDensity();
// radiation length
radl = material->GetRadlen();
- radl /= TG3Units::Length();
+ radl /= TG4G3Units::Length();
absl = 0.; // this parameter is not defined in Geant4
return nofElements;
//
// ---
- CheckStep();
-
- G4ThreeVector theGlobalPoint(xm[0],xm[1],xm[2]);
+ G4AffineTransform affineTransform;
- //const G4NavigationHistory* history
- // = fStep->GetPreStepPoint()->GetTouchable()->GetHistory();
- G4AffineTransform affineTransform
- = fStep->GetPreStepPoint()->GetTouchable()->GetHistory()
- ->GetTopTransform();
+ if (fStepStatus == kVertex) {
+ G4Navigator* navigator =
+ G4TransportationManager::GetTransportationManager()->
+ GetNavigatorForTracking();
+
+ affineTransform = navigator->GetGlobalToLocalTransform();
+ }
+ else {
+
+#ifdef TGEANT4_DEBUG
+ CheckStep("Gmtod");
+#endif
+
+ affineTransform
+ = fStep->GetPreStepPoint()->GetTouchable()->GetHistory()
+ ->GetTopTransform();
+ }
+ G4ThreeVector theGlobalPoint(xm[0],xm[1],xm[2]);
G4ThreeVector theLocalPoint;
if(iflag == 1)
theLocalPoint = affineTransform.TransformPoint(theGlobalPoint);
xd[0] = theLocalPoint.x();
xd[1] = theLocalPoint.y();
xd[2] = theLocalPoint.z();
-
- /*
- // does not work ???
- G4ThreeVector direction(0,0,0);
- G4bool RelativeSearch = true;
- G4Navigator* theNavigator =
- G4TransportationManager::GetTransportationManager()->
- GetNavigatorForTracking();
-
- G4VPhysicalVolume* pPhysVol;
- pPhysVol
- = LocateGlobalPointAndSetup(theGlobalPoint, &direction, RelativeSearch);
- //LocateGlobalPointWithinVolume(theGlobalPoint);
-
- G4AffineTransform at
- = theNavigator->GetGlobalToLocalTransform();
- */
}
void TG4StepManager::Gdtom(Float_t* xd, Float_t* xm, Int_t iflag)
//
// ---
- CheckStep();
+ G4AffineTransform affineTransform;
- // check this
-
- G4ThreeVector theLocalPoint(xd[0],xd[1],xd[2]);
+ if (fStepStatus == kVertex) {
+ G4Navigator* navigator =
+ G4TransportationManager::GetTransportationManager()->
+ GetNavigatorForTracking();
+
+ affineTransform = navigator->GetLocalToGlobalTransform();
+ }
+ else {
- G4AffineTransform affineTransform
- = fStep->GetPreStepPoint()->GetTouchable()->GetHistory()
+#ifdef TGEANT4_DEBUG
+ CheckStep("Gdtom");
+#endif
+
+ // check this
+
+ affineTransform
+ = fStep->GetPreStepPoint()->GetTouchable()->GetHistory()
->GetTopTransform().Inverse();
+ }
+ G4ThreeVector theLocalPoint(xd[0],xd[1],xd[2]);
G4ThreeVector theGlobalPoint;
if(iflag == 1)
theGlobalPoint = affineTransform.TransformPoint(theLocalPoint);
// by User Limits.
// ---
- // check this
- G4LogicalVolume* curLogVolume
+ G4LogicalVolume* curLogVolume
= GetCurrentPhysicalVolume()->GetLogicalVolume();
+
+ // check this
G4UserLimits* userLimits
= curLogVolume->GetUserLimits();
else {
const G4Track& trackRef = *(fTrack);
maxStep = userLimits->GetMaxAllowedStep(trackRef);
- maxStep /= TG3Units::Length();
+ maxStep /= TG4G3Units::Length();
return maxStep;
}
}
// (position of the PostStepPoint).
// ---
+#ifdef TGEANT4_DEBUG
CheckTrack();
-
+#endif
+
// 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);
}
= GetCurrentPhysicalVolume()->GetLogicalVolume()->GetMaterial();
// medium index
- TG4GeometryManager* pGeometryManager = TG4GeometryManager::Instance();
- return pGeometryManager->GetMediumId(curMaterial);
+ TG4GeometryServices* geometryServices = TG4GeometryServices::Instance();
+ return geometryServices->GetMediumId(curMaterial);
}
void TG4StepManager::TrackMomentum(TLorentzVector& momentum) const
// Current particle "momentum" (px, py, pz, Etot).
// ---
+#ifdef TGEANT4_DEBUG
CheckTrack();
+#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);
}
// and the local time since the current track is created.
// ---
+#ifdef TGEANT4_DEBUG
CheckTrack();
+#endif
// 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);
}
// to do: change Ekin -> Etot
// ---
+#ifdef TGEANT4_DEBUG
CheckTrack();
+#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);
}
G4double length;
if (fStepStatus == kNormalStep) {
- CheckStep();
+
+#ifdef TGEANT4_DEBUG
+ CheckStep("TrackStep");
+#endif
+
length = fStep->GetStepLength();
- length /= TG3Units::Length();
+ length /= TG4G3Units::Length();
}
else
length = 0;
// Returns the length of the current track from its origin.
// ---
+#ifdef TGEANT4_DEBUG
CheckTrack();
+#endif
G4double length = fTrack->GetTrackLength();
- length /= TG3Units::Length();
+ length /= TG4G3Units::Length();
return length;
}
// the proper time of the dynamical particle of the current track.
// ---
+#ifdef TGEANT4_DEBUG
CheckTrack();
+#endif
G4double time = fTrack->GetLocalTime();
- time /= TG3Units::Time();
+ time /= TG4G3Units::Time();
return time;
}
G4double energyDeposit;
if (fStepStatus == kNormalStep) {
- CheckStep();
+
+#ifdef TGEANT4_DEBUG
+ CheckStep("Edep");
+#endif
+
energyDeposit = fStep->GetTotalEnergyDeposit();
- energyDeposit /= TG3Units::Energy();
+ energyDeposit /= TG4G3Units::Energy();
}
else
energyDeposit = 0;
// Returns the current particle PDG encoding.
// ---
+#ifdef TGEANT4_DEBUG
CheckTrack();
+#endif
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;
}
// Returns the current particle charge.
// ---
+#ifdef TGEANT4_DEBUG
CheckTrack();
+#endif
+
G4double charge
= fTrack->GetDynamicParticle()->GetDefinition()
->GetPDGCharge();
- charge /= TG3Units::Charge();
+ charge /= TG4G3Units::Charge();
return charge;
}
// Returns current particle rest mass.
// ---
+#ifdef TGEANT4_DEBUG
CheckTrack();
+#endif
G4double mass
= fTrack->GetDynamicParticle()->GetDefinition()
->GetPDGMass();
- mass /= TG3Units::Mass();
+ mass /= TG4G3Units::Mass();
return mass;
}
// Returns total energy of the current particle.
// ---
+#ifdef TGEANT4_DEBUG
CheckTrack();
+#endif
G4double energy
= fTrack->GetDynamicParticle()->GetTotalEnergy();
- energy /= TG3Units::Energy();
+ energy /= TG4G3Units::Energy();
return energy;
}
// ---
if (fStepStatus == kNormalStep) {
- CheckStep();
-
+
+#ifdef TGEANT4_DEBUG
+ CheckStep("IsTrackExiting");
+#endif
+
if (fStep->GetPostStepPoint()->GetStepStatus() == fGeomBoundary)
return true;
}
if (fStepStatus == kVertex) return false;
- CheckStep();
+#ifdef TGEANT4_DEBUG
+ CheckStep("IsTrackCut");
+#endif
// check
G4StepStatus status
// // to the next event.
// ---
+#ifdef TGEANT4_DEBUG
CheckTrack();
+#endif
// check
G4TrackStatus status
// or has been killed, suspended or postponed to next event.
// ---
+#ifdef TGEANT4_DEBUG
CheckTrack();
+#endif
// check
G4TrackStatus status
// Returns true if particle continues tracking.
// ---
+#ifdef TGEANT4_DEBUG
CheckTrack();
+#endif
G4TrackStatus status
= fTrack->GetTrackStatus();
// in the current step.
// ---
+#ifdef TGEANT4_DEBUG
CheckSteppingManager();
+#endif
G4int nofSecondaries = 0;
nofSecondaries += fSteppingManager->GetfN2ndariesAtRestDoIt();
// !! Check if indexing of secondaries is same !!
// ---
+#ifdef TGEANT4_DEBUG
CheckSteppingManager();
+#endif
G4int nofSecondaries = NSecondaries();
G4TrackVector* secondaryTracks = fSteppingManager->GetSecondary();
// position & time
G4ThreeVector positionVector = track->GetPosition();
- positionVector *= 1./(TG3Units::Length());
+ positionVector *= 1./(TG4G3Units::Length());
G4double time = track->GetLocalTime();
- time /= TG3Units::Time();
+ time /= TG4G3Units::Time();
SetTLorentzVector(positionVector, time, position);
// momentum & energy
G4ThreeVector momentumVector = track->GetMomentum();
G4double energy = track->GetDynamicParticle()->GetTotalEnergy();
- energy /= TG3Units::Energy();
+ energy /= TG4G3Units::Energy();
SetTLorentzVector(momentumVector, energy, momentum);
}
else {
}
}
-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.
// ---
- CheckStep();
+ G4int nofSecondaries = NSecondaries();
+ if (fStepStatus == kVertex || !nofSecondaries) return kPNoProcess;
+
+#ifdef TGEANT4_DEBUG
+ CheckStep("ProdProcess");
+#endif
- const G4VProcess* curProcess
- = fStep->GetPostStepPoint()->GetProcessDefinedStep();
+ G4TrackVector* secondaryTracks = fSteppingManager->GetSecondary();
+
+ // should never happen
+ if (!secondaryTracks) {
+ TG4Globals::Exception(
+ "TG4StepManager::ProdProcess(): secondary tracks vector is empty.");
- G4String g4Name = curProcess->GetProcessName();
- return g4Name;
+ return kPNoProcess;
+ }
+
+ if (isec < 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
+
+ // 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;
+ }
+ else {
+ TG4Globals::Exception(
+ "TG4StepManager::GetSecondary(): wrong secondary track index.");
+
+ return kPNoProcess;
+ }
+}
+
+
+Int_t TG4StepManager::StepProcesses(TArrayI &proc) 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) {
+ G4cout << "kVertex" << G4endl;
+ G4int nofProcesses = 1;
+ proc.Set(nofProcesses);
+ proc[0] = kPNull;
+ return nofProcesses;
+ }
+
+#ifdef TGEANT4_DEBUG
+ CheckSteppingManager();
+ CheckStep("StepProcesses");
+#endif
+
+ // along step processes
+ G4ProcessManager* processManager
+ = fStep->GetTrack()->GetDefinition()->GetProcessManager();
+ G4ProcessVector* alongStepProcessVector
+ = processManager->GetAlongStepProcessVector();
+ G4int nofProcesses = alongStepProcessVector->entries();
+
+ // process defined step
+ const G4VProcess* kpLastProcess
+ = fStep->GetPostStepPoint()->GetProcessDefinedStep();
+
+ // fill the array of processes
+ proc.Set(nofProcesses);
+ TG4PhysicsManager* physicsManager = TG4PhysicsManager::Instance();
+ G4int i;
+ for (i=0; i<nofProcesses-1; i++) {
+ G4VProcess* g4Process = (*alongStepProcessVector)[i];
+ // do not fill transportation along step process
+ if (g4Process->GetProcessName() != "Transportation") {
+ physicsManager->GetMCProcess(g4Process);
+ proc[i] = physicsManager->GetMCProcess(g4Process);
+ }
+ }
+ proc[nofProcesses-1] = physicsManager->GetMCProcess(kpLastProcess);
+
+ return nofProcesses;
}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff