Application/src/ElectronLimiter.cc

0001 //
0002 // V.Ivanchenko 2013/10/19
0003 // step limiter and killer for e+,e- and other charged particles
0004 //
0005 #include "SimG4Core/Application/interface/ElectronLimiter.h"
0006 #include "SimG4Core/PhysicsLists/interface/CMSTrackingCutModel.h"
0007 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0008
0009 #include "G4ParticleDefinition.hh"
0010 #include "G4VEnergyLossProcess.hh"
0011 #include "G4LossTableManager.hh"
0012 #include "G4DummyModel.hh"
0013 #include "G4Step.hh"
0014 #include "G4Track.hh"
0015 #include "G4Region.hh"
0016 #include <CLHEP/Units/SystemOfUnits.h>
0017 #include "G4TransportationProcessType.hh"
0018
0019 ElectronLimiter::ElectronLimiter(const edm::ParameterSet &p, const G4ParticleDefinition *part)
0020     : G4VEmProcess("eLimiter", fGeneral),
0021       ionisation_(nullptr),
0022       particle_(part),
0023       nRegions_(0),
0024       rangeCheckFlag_(false),
0025       fieldCheckFlag_(false),
0026       killTrack_(false) {
0027   // set Process Sub Type
0028   SetProcessSubType(static_cast<int>(STEP_LIMITER));
0029   minStepLimit_ = p.getParameter<double>("MinStepLimit") * CLHEP::mm;
0030   trcut_ = new CMSTrackingCutModel(part);
0031 }
0032
0033 ElectronLimiter::~ElectronLimiter() { delete trcut_; }
0034
0035 void ElectronLimiter::InitialiseProcess(const G4ParticleDefinition *) {
0036   G4LossTableManager *lManager = G4LossTableManager::Instance();
0037   if (rangeCheckFlag_) {
0038     ionisation_ = lManager->GetEnergyLossProcess(particle_);
0039     if (!ionisation_) {
0040       rangeCheckFlag_ = false;
0041     }
0042   }
0043   AddEmModel(0, new G4DummyModel());
0044
0045   if (lManager->IsMaster()) {
0046     edm::LogVerbatim("SimG4CoreApplication")
0047         << "ElectronLimiter::BuildPhysicsTable for " << particle_->GetParticleName() << " ioni: " << ionisation_
0048         << " rangeCheckFlag: " << rangeCheckFlag_ << " fieldCheckFlag: " << fieldCheckFlag_ << " " << nRegions_
0049         << " regions for tracking cuts\n";
0050   }
0051 }
0052
0053 G4double ElectronLimiter::PostStepGetPhysicalInteractionLength(const G4Track &aTrack,
0054                                                                G4double previousLimit,
0055                                                                G4ForceCondition *condition) {
0056   *condition = NotForced;
0057   G4double limit = DBL_MAX;
0058   killTrack_ = false;
0059
0060   G4double kinEnergy = aTrack.GetKineticEnergy();
0061   if (0 < nRegions_) {
0062     if (regions_[0] == nullptr) {
0063       if (kinEnergy < energyLimits_[0]) {
0064         killTrack_ = true;
0065         trcut_->InitialiseForStep(factors_[0], rms_[0]);
0066         limit = 0.0;
0067       }
0068     } else {
0069       const G4Region *reg = aTrack.GetVolume()->GetLogicalVolume()->GetRegion();
0070       for (G4int i = 0; i < nRegions_; ++i) {
0071         if (reg == regions_[i] && kinEnergy < energyLimits_[i]) {
0072           killTrack_ = true;
0073           trcut_->InitialiseForStep(factors_[i], rms_[i]);
0074           limit = 0.0;
0075           break;
0076         }
0077       }
0078     }
0079   }
0080   if (!killTrack_ && rangeCheckFlag_) {
0081     G4double safety = aTrack.GetStep()->GetPreStepPoint()->GetSafety();
0082     if (safety > std::min(minStepLimit_, previousLimit)) {
0083       G4double range = ionisation_->GetRange(kinEnergy, aTrack.GetMaterialCutsCouple());
0084       if (safety >= range) {
0085         killTrack_ = true;
0086         limit = 0.0;
0087       }
0088     }
0089   }
0090   if (!killTrack_ && fieldCheckFlag_) {
0091     limit = minStepLimit_;
0092   }
0093   return limit;
0094 }
0095
0096 G4VParticleChange *ElectronLimiter::PostStepDoIt(const G4Track &aTrack, const G4Step &) {
0097   fParticleChange.Initialize(aTrack);
0098   if (killTrack_) {
0099     fParticleChange.ProposeTrackStatus(fStopAndKill);
0100     G4double edep = trcut_->SampleEnergyDepositEcal(aTrack.GetKineticEnergy());
0101     fParticleChange.ProposeLocalEnergyDeposit(edep);
0102     fParticleChange.SetProposedKineticEnergy(0.0);
0103   }
0104   return &fParticleChange;
0105 }
0106
0107 G4bool ElectronLimiter::IsApplicable(const G4ParticleDefinition &) { return true; }
0108
0109 void ElectronLimiter::StartTracking(G4Track *) {}
0110
0111 void ElectronLimiter::SetTrackingCutPerRegion(std::vector<const G4Region *> &reg,
0112                                               std::vector<G4double> &cut,
0113                                               std::vector<G4double> &fac,
0114                                               std::vector<G4double> &rms) {
0115   nRegions_ = reg.size();
0116   if (nRegions_ > 0) {
0117     regions_ = reg;
0118     energyLimits_ = cut;
0119     factors_ = fac;
0120     rms_ = rms;
0121   }
0122 }