]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - TGeant4/TG4StepManager.cxx
IsEmpty method implemented
[u/mrichter/AliRoot.git] / TGeant4 / TG4StepManager.cxx
index 15132822372313c5e1fe86031fda23c520228d9f..5e0e34339796035b8b6117149ce3baf51d7854ef 100644 (file)
@@ -1,29 +1,37 @@
 // $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 <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;
 
+//_____________________________________________________________________________
 TG4StepManager::TG4StepManager() 
   : fTrack(0),
     fStep(0),
@@ -39,18 +47,21 @@ TG4StepManager::TG4StepManager()
   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
@@ -64,6 +75,7 @@ TG4StepManager& TG4StepManager::operator=(const TG4StepManager& right)
           
 // private methods
 
+//_____________________________________________________________________________
 void TG4StepManager::CheckTrack() const
 {
 // Gives exception in case the track is not defined.
@@ -74,16 +86,21 @@ void TG4StepManager::CheckTrack() const
 }     
 
 
-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);
+  }
 }     
 
 
+//_____________________________________________________________________________
 void TG4StepManager::CheckSteppingManager() const
 {
 // Gives exception in case the step is not defined.
@@ -94,6 +111,7 @@ void TG4StepManager::CheckSteppingManager() const
 }     
 
 
+//_____________________________________________________________________________
 void TG4StepManager::SetTLorentzVector(G4ThreeVector xyz, G4double t, 
                                        TLorentzVector& lv) const                                      
 {
@@ -106,6 +124,7 @@ void TG4StepManager::SetTLorentzVector(G4ThreeVector xyz, G4double t,
    lv[3] = t;
 }                                     
 
+//_____________________________________________________________________________
 G4VPhysicalVolume* TG4StepManager::GetCurrentOffPhysicalVolume(G4int off) const 
 {
 // Returns the physical volume of the off-th mother's
@@ -114,7 +133,8 @@ G4VPhysicalVolume* TG4StepManager::GetCurrentOffPhysicalVolume(G4int off) const
  
   G4VPhysicalVolume* physVolume = GetCurrentPhysicalVolume(); 
 
-  G4VPhysicalVolume* mother = physVolume;
+  G4VPhysicalVolume* mother = physVolume; 
+
   Int_t level = off;
   while (level > 0) { 
     if (mother) mother = mother->GetMother();
@@ -132,38 +152,9 @@ G4VPhysicalVolume* TG4StepManager::GetCurrentOffPhysicalVolume(G4int off) const
   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.
@@ -178,19 +169,24 @@ void TG4StepManager::StopTrack()
 //                       // secondaries.
 // ---
 
+#ifdef TGEANT4_DEBUG
   CheckTrack();
+#endif  
   
   fTrack->SetTrackStatus(fStopAndKill);
   // fTrack->SetTrackStatus(fStopButAlive);
   // fTrack->SetTrackStatus(fKillTrackAndSecondaries);
 }
 
+//_____________________________________________________________________________
 void TG4StepManager::StopEvent()
 {
 // Aborts the current event processing.
 // ---
 
+#ifdef TGEANT4_DEBUG
   CheckTrack();
+#endif  
   
   fTrack->SetTrackStatus(fKillTrackAndSecondaries);
           //StopTrack();   // cannot be used as it keeps secondaries
@@ -198,38 +194,30 @@ void TG4StepManager::StopEvent()
   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.
@@ -239,6 +227,7 @@ void TG4StepManager::SetMaxNStep(Int_t maxNofSteps)
     "TG4StepManager::SetMaxNStep(..) is not yet implemented.");
 }
 
+//_____________________________________________________________________________
 void TG4StepManager::SetUserDecay(Int_t pdg)
 {
 // Not yet implemented.
@@ -248,6 +237,7 @@ void TG4StepManager::SetUserDecay(Int_t pdg)
     "TG4StepManager::SetUserDecay(..) is not yet implemented.");
 }
 
+//_____________________________________________________________________________
 G4VPhysicalVolume* TG4StepManager::GetCurrentPhysicalVolume() const 
 {
 // Returns the current physical volume.
@@ -257,15 +247,27 @@ G4VPhysicalVolume* TG4StepManager::GetCurrentPhysicalVolume() const
 
   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()->
@@ -277,6 +279,7 @@ G4VPhysicalVolume* TG4StepManager::GetCurrentPhysicalVolume() const
   return physVolume;
 }     
 
+//_____________________________________________________________________________
 Int_t TG4StepManager::CurrentVolID(Int_t& copyNo) const
 {
 // Returns the current sensitive detector ID
@@ -287,9 +290,11 @@ Int_t TG4StepManager::CurrentVolID(Int_t& copyNo) const
   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
@@ -304,7 +309,8 @@ Int_t TG4StepManager::CurrentVolOffID(Int_t off, Int_t&  copyNo) const
     copyNo = mother->GetCopyNo() + 1;
 
     // sensitive detector ID
-    return GetVolumeID(mother);
+    TG4SDServices* sdServices = TG4SDServices::Instance();
+    return sdServices->GetVolumeID(mother->GetLogicalVolume());
   }
   else {
     copyNo = 0;
@@ -312,6 +318,7 @@ Int_t TG4StepManager::CurrentVolOffID(Int_t off, Int_t&  copyNo) const
   }  
 }
 
+//_____________________________________________________________________________
 const char* TG4StepManager::CurrentVolName() const
 {
 // Returns the current physical volume name.
@@ -320,6 +327,7 @@ const char* TG4StepManager::CurrentVolName() const
   return GetCurrentPhysicalVolume()->GetName();
 }
 
+//_____________________________________________________________________________
 const char* TG4StepManager::CurrentVolOffName(Int_t off) const
 { 
 // Returns the off-th mother's physical volume name.
@@ -335,6 +343,7 @@ const char* TG4StepManager::CurrentVolOffName(Int_t off) const
     return 0;
 }
 
+//_____________________________________________________________________________
 Int_t TG4StepManager::CurrentMaterial(Float_t &a, Float_t &z, Float_t &dens, 
                           Float_t &radl, Float_t &absl) const
 {
@@ -347,32 +356,32 @@ Int_t TG4StepManager::CurrentMaterial(Float_t &a, Float_t &z, Float_t &dens,
   G4Material* material 
     = physVolume->GetLogicalVolume()->GetMaterial();
 
-  // this if may be redundant - check
-  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
@@ -393,31 +402,48 @@ void TG4StepManager::Gmtod(Float_t* xm, Float_t* xd, Int_t iflag)
 //
 // ---
 
-  CheckStep();
-  G4ThreeVector theGlobalPoint(xm[0],xm[1],xm[2]); 
+#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) {
+    G4Navigator* navigator =
+      G4TransportationManager::GetTransportationManager()->
+        GetNavigatorForTracking();
+       
+    affineTransform = navigator->GetGlobalToLocalTransform();
+  }
+  else {
 
-  //const G4NavigationHistory* history
-  //   =  fStep->GetPreStepPoint()->GetTouchable()->GetHistory();
-  G4AffineTransform affineTransform
-    = fStep->GetPreStepPoint()->GetTouchable()->GetHistory()
-      ->GetTopTransform();
+#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);
-  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
@@ -437,16 +463,29 @@ 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 {
+
+#ifdef TGEANT4_DEBUG
+    CheckStep("Gdtom");
+#endif
 
-  G4AffineTransform affineTransform
-    = fStep->GetPreStepPoint()->GetTouchable()->GetHistory()
+    // 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);
@@ -461,15 +500,17 @@ void TG4StepManager::Gdtom(Float_t* xd, Float_t* xm, Int_t iflag)
   xm[2] = theGlobalPoint.z();
 } 
  
+//_____________________________________________________________________________
 Float_t TG4StepManager::MaxStep() const
 {   
 // Returns maximum step allowed in the current logical volume
 // by User Limits.
 // ---
 
-  // check this
-  G4LogicalVolume* curLogVolume 
+  G4LogicalVolume* curLogVolume
     = GetCurrentPhysicalVolume()->GetLogicalVolume();
+
+  // check this
   G4UserLimits* userLimits 
     = curLogVolume->GetUserLimits();
 
@@ -483,11 +524,12 @@ Float_t TG4StepManager::MaxStep() const
   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.
@@ -498,6 +540,7 @@ Int_t TG4StepManager::GetMaxNStep() const
   return 0; 
 }
 
+//_____________________________________________________________________________
 void TG4StepManager::TrackPosition(TLorentzVector& position) const
 { 
 // Current particle position (in the world reference frame)
@@ -505,87 +548,100 @@ void TG4StepManager::TrackPosition(TLorentzVector& position) const
 // (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);
 }
 
+//_____________________________________________________________________________
 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).
 // ---
 
+#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);
 }
 
+//_____________________________________________________________________________
 void TG4StepManager::TrackVertexPosition(TLorentzVector& position) const
 { 
 // The vertex particle position (in the world reference frame)
 // 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);
 }
 
+//_____________________________________________________________________________
 void TG4StepManager::TrackVertexMomentum(TLorentzVector& momentum) const
 {  
 // The vertex particle "momentum" (px, py, pz, Ekin)
 // 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);
 }
 
+//_____________________________________________________________________________
 Float_t TG4StepManager::TrackStep() const
 {   
 // Returns the current step length.
@@ -593,9 +649,13 @@ Float_t TG4StepManager::TrackStep() const
 
   G4double length;
   if (fStepStatus == kNormalStep) {
-    CheckStep();    
+
+#ifdef TGEANT4_DEBUG
+    CheckStep("TrackStep");    
+#endif
+
     length = fStep->GetStepLength();
-    length /= TG3Units::Length();
+    length /= TG4G3Units::Length();
   }  
   else 
     length = 0;
@@ -603,18 +663,22 @@ Float_t TG4StepManager::TrackStep() const
   return length;
 }
 
+//_____________________________________________________________________________
 Float_t TG4StepManager::TrackLength() const
 {
 // 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;
 }
 
+//_____________________________________________________________________________
 Float_t TG4StepManager::TrackTime() const
 {
 // Returns the local time since the current track is created.
@@ -623,13 +687,16 @@ Float_t TG4StepManager::TrackTime() const
 // 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;
 }
 
+//_____________________________________________________________________________
 Float_t TG4StepManager::Edep() const
 {   
 // Returns total energy deposit in this step.
@@ -637,9 +704,13 @@ Float_t TG4StepManager::Edep() const
 
   G4double energyDeposit;
   if (fStepStatus == kNormalStep) {
-    CheckStep();
+
+#ifdef TGEANT4_DEBUG
+    CheckStep("Edep");
+#endif
+
     energyDeposit = fStep->GetTotalEnergyDeposit();
-    energyDeposit /= TG3Units::Energy();
+    energyDeposit /= TG4G3Units::Energy();
   }
   else   
     energyDeposit = 0;
@@ -647,65 +718,79 @@ Float_t TG4StepManager::Edep() const
   return energyDeposit;
 }
 
+//_____________________________________________________________________________
 Int_t TG4StepManager::TrackPid() const
 {   
 // 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;
 }
 
+//_____________________________________________________________________________
 Float_t TG4StepManager::TrackCharge() const
 {   
 // Returns the current particle charge.
 // ---
 
+#ifdef TGEANT4_DEBUG
   CheckTrack();
+#endif
+
   G4double charge
     = fTrack->GetDynamicParticle()->GetDefinition()
       ->GetPDGCharge();
-  charge /= TG3Units::Charge();        
+  charge /= TG4G3Units::Charge();      
   return charge;
 }
 
+//_____________________________________________________________________________
 Float_t TG4StepManager::TrackMass() const
 {   
 // Returns current particle rest mass.
 // ---
 
+#ifdef TGEANT4_DEBUG
   CheckTrack();
+#endif
 
   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.
 // ---
 
+#ifdef TGEANT4_DEBUG
   CheckTrack();
+#endif
 
   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
@@ -720,6 +805,7 @@ Bool_t TG4StepManager::IsTrackInside() const
   return false;    
 }
 
+//_____________________________________________________________________________
 Bool_t TG4StepManager::IsTrackEntering() const
 {   
 // Returns true if particle cross a geometrical boundary
@@ -734,14 +820,18 @@ Bool_t TG4StepManager::IsTrackEntering() const
   return false;  
 }
 
+//_____________________________________________________________________________
 Bool_t TG4StepManager::IsTrackExiting() const
 {   
 // Returns true if particle cross a geometrical boundary.
 // ---
 
   if (fStepStatus == kNormalStep) {
-    CheckStep();
-    
+
+#ifdef TGEANT4_DEBUG
+    CheckStep("IsTrackExiting");
+#endif    
+
     if (fStep->GetPostStepPoint()->GetStepStatus() == fGeomBoundary) 
        return true;  
   }
@@ -749,6 +839,7 @@ Bool_t TG4StepManager::IsTrackExiting() const
   return false;  
 }
 
+//_____________________________________________________________________________
 Bool_t TG4StepManager::IsTrackOut() const
 {   
 // Returns true if particle cross the world boundary
@@ -757,7 +848,9 @@ Bool_t TG4StepManager::IsTrackOut() const
 
   if (fStepStatus == kVertex) return false;
 
-  CheckStep();
+#ifdef TGEANT4_DEBUG
+  CheckStep("IsTrackCut");
+#endif
 
   // check
   G4StepStatus status
@@ -768,6 +861,7 @@ Bool_t TG4StepManager::IsTrackOut() const
     return false;
 }
 
+//_____________________________________________________________________________
 Bool_t TG4StepManager::IsTrackStop() const
 {   
 // Returns true if particle has stopped 
@@ -787,7 +881,9 @@ Bool_t TG4StepManager::IsTrackStop() const
 //                       // to the next event.
 // ---
 
+#ifdef TGEANT4_DEBUG
   CheckTrack();
+#endif
 
   // check
   G4TrackStatus status
@@ -802,6 +898,7 @@ Bool_t TG4StepManager::IsTrackStop() const
     return false; 
 }
 
+//_____________________________________________________________________________
 Bool_t TG4StepManager::IsTrackDisappeared() const
 { 
 // Returns true if particle has disappeared 
@@ -809,7 +906,9 @@ Bool_t TG4StepManager::IsTrackDisappeared() const
 // or has been killed, suspended or postponed to next event.
 // ---
 
+#ifdef TGEANT4_DEBUG
   CheckTrack();
+#endif
 
   // check
   G4TrackStatus status
@@ -824,12 +923,15 @@ Bool_t TG4StepManager::IsTrackDisappeared() const
     return false;
 }
 
+//_____________________________________________________________________________
 Bool_t TG4StepManager::IsTrackAlive() const
 {   
 // Returns true if particle continues tracking.
 // ---
 
+#ifdef TGEANT4_DEBUG
   CheckTrack();
+#endif
 
   G4TrackStatus status
      = fTrack->GetTrackStatus();
@@ -839,6 +941,7 @@ Bool_t TG4StepManager::IsTrackAlive() const
     return false; 
 }
 
+//_____________________________________________________________________________
 Bool_t TG4StepManager::IsNewTrack() const
 {
 // Returns true when track performs the first step.
@@ -850,13 +953,16 @@ Bool_t TG4StepManager::IsNewTrack() const
     return false;
 }
 
+//_____________________________________________________________________________
 Int_t TG4StepManager::NSecondaries() const
 {
 // Returns the number of secondary particles generated 
 // in the current step.
 // ---
 
+#ifdef TGEANT4_DEBUG
   CheckSteppingManager();
+#endif
 
   G4int nofSecondaries = 0;
   nofSecondaries += fSteppingManager->GetfN2ndariesAtRestDoIt();
@@ -866,6 +972,7 @@ Int_t TG4StepManager::NSecondaries() const
   return nofSecondaries;
 }
 
+//_____________________________________________________________________________
 void TG4StepManager::GetSecondary(Int_t index, Int_t& particleId, 
                           TLorentzVector& position, TLorentzVector& momentum)
 {
@@ -874,61 +981,157 @@ void TG4StepManager::GetSecondary(Int_t index, Int_t& particleId,
 // !! Check if indexing of secondaries is same !!
 // ---
 
+#ifdef TGEANT4_DEBUG
   CheckSteppingManager();
+#endif
 
   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.
 // ---
 
-  CheckStep();
+  G4int nofSecondaries = NSecondaries();
+  if (fStepStatus == kVertex || !nofSecondaries) return kPNoProcess;
+
+#ifdef TGEANT4_DEBUG
+  CheckStep("ProdProcess");
+#endif
+
+  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.");
+
+    return kPNoProcess;  
+  }
+#endif
 
-  const G4VProcess* curProcess 
-    = fStep->GetPostStepPoint()->GetProcessDefinedStep(); 
+  // 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
 
-  G4String g4Name = curProcess->GetProcessName(); 
-  return g4Name;
+  // 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;  
 }