TGeant4/TG4SpecialCuts.cxx

   1 // $Id$
   2 // Category: physics
   3 //
   4 // See the class description in the header file.
   5
   6 #include "TG4SpecialCuts.h"
   7 #include "TG4CutVector.h"
   8 #include "TG4Limits.h"
   9
  10 #include <G4UserLimits.hh>
  11
  12 #include <G4EnergyLossTables.hh>
  13
  14 TG4SpecialCuts::TG4SpecialCuts(TG3ParticleWSP particle, TG4CutVector* cutVector,
  15                                const G4String& processName)
  16   : G4UserSpecialCuts(processName),
  17     fCutVector(cutVector)
  18 {
  19 //
  20   switch (particle) {
  21     case kGamma:
  22       fPtrMinEkineInCutVector = &TG4CutVector::GetMinEkineForGamma;
  23       fPtrMinEkineInLimits = &TG4Limits::GetMinEkineForGamma;
  24       break;
  25     case kElectron: case kEplus:
  26       fPtrMinEkineInCutVector = &TG4CutVector::GetMinEkineForElectron;
  27       fPtrMinEkineInLimits = &TG4Limits::GetMinEkineForElectron;
  28       break;
  29     case kChargedHadron:
  30       fPtrMinEkineInCutVector = &TG4CutVector::GetMinEkineForHadron;
  31       fPtrMinEkineInLimits = &TG4Limits::GetMinEkineForHadron;
  32       break;
  33     case kNeutralHadron:
  34       fPtrMinEkineInCutVector = &TG4CutVector::GetMinEkineForNeutralHadron;
  35       fPtrMinEkineInLimits = &TG4Limits::GetMinEkineForNeutralHadron;
  36       break;
  37     case kMuon:
  38       fPtrMinEkineInCutVector = &TG4CutVector::GetMinEkineForMuon;
  39       fPtrMinEkineInLimits = &TG4Limits::GetMinEkineForMuon;
  40       break;
  41     case kAny:
  42       fPtrMinEkineInCutVector = &TG4CutVector::GetMinEkineForOther;
  43       fPtrMinEkineInLimits = &TG4Limits::GetMinEkineForOther;
  44       break;
  45     case kNofParticlesWSP:
  46       TG4Globals::Exception("TG4SpecialCuts: Wrong particle specification.");
  47       break;
  48   }
  49 }
  50
  51 TG4SpecialCuts::TG4SpecialCuts() {
  52 //
  53 }
  54
  55 TG4SpecialCuts::TG4SpecialCuts(const TG4SpecialCuts& right) {
  56 //
  57   TG4Globals::Exception(
  58     "TG4SpecialCuts is protected from copying.");
  59 }
  60
  61 TG4SpecialCuts::~TG4SpecialCuts() {
  62 //
  63 }
  64
  65 // operators
  66
  67 TG4SpecialCuts& TG4SpecialCuts::operator=(const TG4SpecialCuts& right)
  68 {
  69   // check assignement to self
  70   if (this == &right) return *this;
  71
  72   TG4Globals::Exception(
  73     "TG4SpecialCuts is protected from assigning.");
  74
  75   return *this;
  76 }
  77
  78 // public methods
  79
  80 G4double TG4SpecialCuts::PostStepGetPhysicalInteractionLength(
  81                            const G4Track& track, G4double previousStepSize,
  82                            G4ForceCondition* condition)
  83 {
  84 // Returns the Step-size (actual length) which is allowed
  85 // by this process.
  86 // ---
  87
  88   // set condition
  89   *condition = NotForced;
  90
  91   G4double proposedStep = DBL_MAX;
  92   G4UserLimits* limits
  93     = track.GetVolume()->GetLogicalVolume()->GetUserLimits();
  94   if (limits) {
  95     // max track length
  96     proposedStep
  97       = (limits->GetUserMaxTrackLength(track) - track.GetTrackLength());
  98     if (proposedStep < 0.) return 0.;
  99
 100     // max time limit
 101     G4double beta
 102       = (track.GetDynamicParticle()->GetTotalMomentum()) /
 103         (track.GetTotalEnergy());
 104     G4double time
 105       = (limits->GetUserMaxTime(track) - track.GetGlobalTime());
 106     G4double temp = beta*c_light*time;
 107     if (temp < 0.) return 0.;
 108     if (proposedStep > temp) proposedStep = temp;
 109
 110     // min remaining range
 111     G4ParticleDefinition* particle = track.GetDefinition();
 112     G4double kinEnergy = track.GetKineticEnergy();
 113     G4Material* material = track.GetMaterial();
 114     G4double rangeNow
 115       = G4EnergyLossTables::GetRange(particle, kinEnergy, material);
 116     temp = (rangeNow - limits->GetUserMinRange(track));
 117     if (temp < 0.) return 0.;
 118     if (proposedStep > temp) proposedStep = temp;
 119
 120     // min kinetic energy (from limits)
 121     // the kin energy cut can be applied only in case
 122     // G4EnergyLossTables are defined for the particle
 123     if (G4EnergyLossTables::GetDEDXTable(particle)) {
 124       TG4Limits* tg4Limits = dynamic_cast<TG4Limits*>(limits);
 125       if (!tg4Limits) {
 126         G4String text = "TG4SpecialCuts::PostStepGetPhysicalInteractionLength:\n";
 127         text = text + "    Unknown limits type.";
 128         TG4Globals::Exception(text);
 129       }
 130       G4double minEkine
 131         = (tg4Limits->*fPtrMinEkineInLimits)(track);
 132       G4double minR
 133         = G4EnergyLossTables::GetRange(particle, minEkine, material);
 134       temp = rangeNow - minR;
 135       if (temp < 0.) return 0.;
 136       if (proposedStep > temp) proposedStep = temp;
 137     }
 138
 139   }
 140   else if (fCutVector) {
 141     // min kinetic energy (from cut vector)
 142     // the kin energy cut can be applied only in case
 143     // G4EnergyLossTables are defined for the particle
 144     G4ParticleDefinition* particle = track.GetDefinition();
 145     if (G4EnergyLossTables::GetDEDXTable(particle)) {
 146       G4double kinEnergy = track.GetKineticEnergy();
 147       G4Material* material = track.GetMaterial();
 148       G4double rangeNow
 149         = G4EnergyLossTables::GetRange(particle, kinEnergy, material);
 150       G4double minEkine
 151         = (fCutVector->*fPtrMinEkineInCutVector)(track);
 152       G4double minR
 153         = G4EnergyLossTables::GetRange(particle, minEkine, material);
 154       G4double temp = rangeNow - minR;
 155       if (temp < 0.) return 0.;
 156       if (proposedStep > temp) proposedStep = temp;
 157     }
 158   }
 159   return proposedStep;
 160 }
 161
 162 G4VParticleChange* TG4SpecialCuts::PostStepDoIt(const G4Track& track,
 163                                                 const G4Step& step)
 164 {
 165 // Kills the current particle, if requested by G4UserLimits.
 166 // ---
 167
 168   aParticleChange.Initialize(track);
 169   aParticleChange.SetEnergyChange(0.) ;
 170   aParticleChange.SetLocalEnergyDeposit(track.GetKineticEnergy()) ;
 171   aParticleChange.SetStatusChange(fStopAndKill);
 172   return &aParticleChange;
 173 }