TGeant4/TG4SpecialCuts.cxx

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