TGeant4/TG4VSpecialCuts.cxx

   1 // $Id$
   2 // Category: physics
   3 //
   4 // Author: I. Hrivnacova
   5 //
   6 // Class TG4VSpecialCuts
   7 // ---------------------
   8 // See the class description in the header file.
   9
  10 #include "TG4VSpecialCuts.h"
  11 #include "TG4G3CutVector.h"
  12 #include "TG4GeometryServices.h"
  13 #include "TG4Limits.h"
  14
  15 #include <G4UserLimits.hh>
  16 #include <G4EnergyLossTables.hh>
  17
  18 //_____________________________________________________________________________
  19 TG4VSpecialCuts::TG4VSpecialCuts(const G4String& processName)
  20   : G4VProcess(processName) {
  21 //
  22 }
  23
  24 //_____________________________________________________________________________
  25 TG4VSpecialCuts::TG4VSpecialCuts() {
  26 //
  27 }
  28
  29 //_____________________________________________________________________________
  30 TG4VSpecialCuts::TG4VSpecialCuts(const TG4VSpecialCuts& right) {
  31 //
  32   TG4Globals::Exception(
  33     "TG4VSpecialCuts is protected from copying.");
  34 }
  35
  36 //_____________________________________________________________________________
  37 TG4VSpecialCuts::~TG4VSpecialCuts() {
  38 //
  39 }
  40
  41 // operators
  42
  43 //_____________________________________________________________________________
  44 TG4VSpecialCuts& TG4VSpecialCuts::operator=(const TG4VSpecialCuts& right)
  45 {
  46   // check assignement to self
  47   if (this == &right) return *this;
  48
  49   TG4Globals::Exception(
  50     "TG4VSpecialCuts is protected from assigning.");
  51
  52   return *this;
  53 }
  54
  55 // public methods
  56
  57 //_____________________________________________________________________________
  58 G4double TG4VSpecialCuts::PostStepGetPhysicalInteractionLength(
  59                            const G4Track& track, G4double previousStepSize,
  60                            G4ForceCondition* condition)
  61 {
  62 // Returns the Step-size (actual length) which is allowed
  63 // by this process.
  64 // ---
  65
  66   // set condition
  67   *condition = NotForced;
  68
  69   G4double proposedStep = DBL_MAX;
  70   G4double minStep = (1.0e-9)*m;
  71
  72   // get limits
  73 #ifdef TGEANT4_DEBUG
  74   TG4Limits* limits
  75      = TG4GeometryServices::Instance()
  76          ->GetLimits(track.GetVolume()->GetLogicalVolume()->GetUserLimits());
  77
  78   if (!limits) {
  79     G4String text = "TG4VSpecialCuts::PostStepGetPhysicalInteractionLength:\n";
  80     text = text + "    " + track.GetVolume()->GetLogicalVolume()->GetName();
  81     text = text + " has not limits.";
  82     TG4Globals::Exception(text);
  83   }
  84 #else
  85   TG4Limits* limits
  86     = (TG4Limits*) track.GetVolume()->GetLogicalVolume()->GetUserLimits();
  87 #endif
  88
  89   // max track length
  90   proposedStep
  91     = (limits->GetUserMaxTrackLength(track) - track.GetTrackLength());
  92   if (proposedStep < 0.) return minStep;
  93
  94   // max time limit
  95   G4double beta
  96     = (track.GetDynamicParticle()->GetTotalMomentum()) /
  97       (track.GetTotalEnergy());
  98   G4double time
  99     = (limits->GetUserMaxTime(track) - track.GetGlobalTime());
 100   G4double temp = beta*c_light*time;
 101   if (temp < 0.) return minStep;
 102   if (proposedStep > temp) proposedStep = temp;
 103
 104   G4ParticleDefinition* particle = track.GetDefinition();
 105   if (particle->GetPDGCharge() != 0.) {
 106
 107     // min remaining range
 108     G4double kinEnergy = track.GetKineticEnergy();
 109     G4Material* material = track.GetMaterial();
 110     G4double rangeNow
 111       = G4EnergyLossTables::GetRange(particle, kinEnergy, material);
 112
 113     temp = (rangeNow - limits->GetUserMinRange(track));
 114     if (temp < 0.) return minStep;
 115     if (proposedStep > temp) proposedStep = temp;
 116
 117     // min kinetic energy (from limits)
 118     // the kin energy cut can be applied only in case
 119     // G4EnergyLossTables are defined for the particle
 120     if (G4EnergyLossTables::GetDEDXTable(particle)) {
 121       G4double minEkine = GetMinEkine(*limits, track);
 122       G4double minR
 123         = G4EnergyLossTables::GetRange(particle, minEkine, material);
 124       temp = rangeNow - minR;
 125       if (temp < 0.) return minStep;
 126       if (proposedStep > temp) proposedStep = temp;
 127     }
 128   }
 129   else {
 130     // min kinetic energy (from limits)
 131     // for neutral particles
 132     G4double minEkine = GetMinEkine(*limits, track);
 133     if (track.GetKineticEnergy() <= minEkine) return minStep;
 134   }
 135
 136   return proposedStep;
 137 }
 138
 139 //_____________________________________________________________________________
 140 G4VParticleChange* TG4VSpecialCuts::PostStepDoIt(const G4Track& track,
 141                                                  const G4Step& step)
 142 {
 143 // Kills the current particle, if requested by G4UserLimits.
 144 // ---
 145
 146   aParticleChange.Initialize(track);
 147   aParticleChange.SetEnergyChange(0.) ;
 148   aParticleChange.SetLocalEnergyDeposit(track.GetKineticEnergy()) ;
 149   aParticleChange.SetStatusChange(fStopAndKill);
 150   return &aParticleChange;
 151 }