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