]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - TGeant4/TG4StepManager.cxx
IsEmpty method implemented
[u/mrichter/AliRoot.git] / TGeant4 / TG4StepManager.cxx
index b1236ffd22d9d74181b4112f5e331f74930ce25d..5e0e34339796035b8b6117149ce3baf51d7854ef 100644 (file)
@@ -1,6 +1,10 @@
 // $Id$
 // Category: event
 //
+// Author: I.Hrivnacova
+//
+// Class TG4StepManager
+// --------------------
 // See the class description in the header file.
 
 #include "TG4StepManager.h"
@@ -9,7 +13,6 @@
 #include "TG4ParticlesManager.h"
 #include "TG4PhysicsManager.h"
 #include "TG4VSensitiveDetector.h"
-#include "TG4Limits.h"
 #include "TG4Globals.h"
 #include "TG4G3Units.h"
 
@@ -197,20 +200,17 @@ void TG4StepManager::SetMaxStep(Float_t step)
 // Maximum step allowed in the current logical volume.
 // ---
 
-  G4LogicalVolume* curLogVolume 
-    = GetCurrentPhysicalVolume()->GetLogicalVolume();
   G4UserLimits* userLimits 
-    = curLogVolume->GetUserLimits();
-
-  if (userLimits == 0) {
-    // 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);
+    = GetCurrentPhysicalVolume()->GetLogicalVolume()->GetUserLimits();
+  
+#ifdef TGEANT4_DEBUG
+  if (!userLimits) {
+    G4String text = "TG4StepManager::SetMaxStep:\n";
+    text = text + "   User limits not defined.";
+    TG4Globals::Exception(text);
+    return;
   }
+#endif  
 
   // set max step
   userLimits->SetMaxAllowedStep(step*TG4G3Units::Length()); 
@@ -356,28 +356,29 @@ Int_t TG4StepManager::CurrentMaterial(Float_t &a, Float_t &z, Float_t &dens,
   G4Material* material 
     = physVolume->GetLogicalVolume()->GetMaterial();
 
-  if (material) {
-    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;
-  }
-  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;
 }
 
 //_____________________________________________________________________________
@@ -401,6 +402,14 @@ void TG4StepManager::Gmtod(Float_t* xm, Float_t* xd, Int_t iflag)
 //
 // ---
 
+#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) {
@@ -423,13 +432,11 @@ void TG4StepManager::Gmtod(Float_t* xm, Float_t* xd, Int_t iflag)
 
   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();
@@ -564,13 +571,12 @@ Int_t TG4StepManager::GetMedium() const
 // G3 tracking medium index).
 // --- 
 
-  // current material
-  G4Material* curMaterial
-    = GetCurrentPhysicalVolume()->GetLogicalVolume()->GetMaterial();
+  // current logical volume
+  G4LogicalVolume* curLV = GetCurrentPhysicalVolume()->GetLogicalVolume();
 
   // medium index  
   TG4GeometryServices* geometryServices = TG4GeometryServices::Instance();
-  return geometryServices->GetMediumId(curMaterial);
+  return geometryServices->GetMediumId(curLV);
 }
 
 //_____________________________________________________________________________
@@ -982,43 +988,42 @@ void TG4StepManager::GetSecondary(Int_t index, Int_t& particleId,
   G4int nofSecondaries = NSecondaries();
   G4TrackVector* secondaryTracks = fSteppingManager->GetSecondary();
 
-  if (secondaryTracks){
-    if (index < nofSecondaries) {
-
-      // 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);
-    }
-    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);
 }
 
 //_____________________________________________________________________________
@@ -1037,6 +1042,7 @@ AliMCProcess TG4StepManager::ProdProcess(Int_t isec) const
 
   G4TrackVector* secondaryTracks = fSteppingManager->GetSecondary();
  
+#ifdef TGEANT4_DEBUG
   // should never happen
   if (!secondaryTracks) {
     TG4Globals::Exception(
@@ -1045,37 +1051,36 @@ AliMCProcess TG4StepManager::ProdProcess(Int_t isec) const
     return kPNoProcess;  
   }    
 
-  if (isec < nofSecondaries) {
+  if (isec >= nofSecondaries) {
+    TG4Globals::Exception(
+      "TG4StepManager::GetSecondary(): wrong secondary track index.");
+
+    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 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]; 
+  // the secondary track with specified isec index
+  G4Track* track = (*secondaryTracks)[startIndex + isec]; 
    
-    const G4VProcess* kpProcess = track->GetCreatorProcess(); 
+  const G4VProcess* kpProcess = track->GetCreatorProcess(); 
   
-    AliMCProcess mcProcess 
-     = TG4PhysicsManager::Instance()->GetMCProcess(kpProcess);
+  AliMCProcess mcProcess 
+   = TG4PhysicsManager::Instance()->GetMCProcess(kpProcess);
   
-    // distinguish kPDeltaRay from kPEnergyLoss  
-    if (mcProcess == kPEnergyLoss) mcProcess = kPDeltaRay;
+  // distinguish kPDeltaRay from kPEnergyLoss  
+  if (mcProcess == kPEnergyLoss) mcProcess = kPDeltaRay;
   
-    return mcProcess;
-  }
-  else {
-    TG4Globals::Exception(
-      "TG4StepManager::GetSecondary(): wrong secondary track index.");
-
-    return kPNoProcess;  
-  }
+  return mcProcess;
 }
 
 //_____________________________________________________________________________
-Int_t TG4StepManager::StepProcesses(TArrayI &proc) const
+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.
@@ -1083,10 +1088,9 @@ Int_t TG4StepManager::StepProcesses(TArrayI &proc) const
 // ---
 
  if (fStepStatus == kVertex) {
-   G4cout << "kVertex" << G4endl;
    G4int nofProcesses = 1;
-   proc.Set(nofProcesses);
-   proc[0] = kPNull;
+   processes.Set(nofProcesses);
+   processes[0] = kPNull;
    return nofProcesses;
  }  
    
@@ -1096,29 +1100,38 @@ Int_t TG4StepManager::StepProcesses(TArrayI &proc) const
 #endif
 
   // along step processes
-  G4ProcessManager* processManager
-    = fStep->GetTrack()->GetDefinition()->GetProcessManager();
-  G4ProcessVector* alongStepProcessVector 
-    = processManager->GetAlongStepProcessVector();
-  G4int nofProcesses = alongStepProcessVector->entries();
+  G4ProcessVector* processVector 
+    = fStep->GetTrack()->GetDefinition()->GetProcessManager()
+        ->GetAlongStepProcessVector();
+  G4int nofAlongStep = processVector->entries();
   
   // process defined step
   const G4VProcess* kpLastProcess 
     = fStep->GetPostStepPoint()->GetProcessDefinedStep();
 
-  // fill the array of processes 
-  proc.Set(nofProcesses);
+  // 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 i;  
-  for (i=0; i<nofProcesses-1; i++) {
-    G4VProcess* g4Process = (*alongStepProcessVector)[i];    
+  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") {
-      physicsManager->GetMCProcess(g4Process);   
-      proc[i] = physicsManager->GetMCProcess(g4Process);
-    }  
-  }  
-  proc[nofProcesses-1] = physicsManager->GetMCProcess(kpLastProcess);
+    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 nofProcesses;  
+  return counter;  
 }