101 currentCoupleIndex = 0;
102 currentCouple =
nullptr;
103 currentMaterial = cutMaterial =
nullptr;
104 currentParticle =
nullptr;
105 lambdaParticle =
nullptr;
106 baseParticle =
nullptr;
107 currentLambda =
nullptr;
108 currentModel =
nullptr;
109 currentProcess =
nullptr;
110 curProcess =
nullptr;
116 cutenergy[0] = cutenergy[1] = cutenergy[2] =
DBL_MAX;
117 currentParticleName=
"";
118 currentMaterialName=
"";
125 isApplicable =
false;
133 for (
G4int i=0; i<nLocalMaterials; ++i) {
134 delete localCouples[i];
147 if(couple && UpdateParticle(p, kinEnergy) ) {
148 res = manager->
GetDEDX(p, kinEnergy, couple);
151 if(FindEmModel(p, currentProcessName, kinEnergy)) {
167 G4cout <<
"G4EmCalculator::GetDEDX: E(MeV)= " << kinEnergy/
MeV
168 <<
" DEDX(MeV/mm)= " << res*
mm/
MeV
172 <<
" isIon= " << isIon
188 if(couple && UpdateParticle(p, kinEnergy)) {
191 G4cout <<
" G4EmCalculator::GetRangeFromRestrictedDEDX: E(MeV)= "
193 <<
" range(mm)= " << res/
mm
212 ed <<
"G4EmCalculator::GetCSDARange: CSDA table is not built; "
213 <<
" use UI command: /process/eLoss/CSDARange true";
214 G4Exception(
"G4EmCalculator::GetCSDARange",
"em0077",
220 if(couple && UpdateParticle(p, kinEnergy)) {
223 G4cout <<
" G4EmCalculator::GetCSDARange: E(MeV)= " << kinEnergy/
MeV
224 <<
" range(mm)= " << res/
mm
258 if(couple && UpdateParticle(p, 1.0*
GeV)) {
259 res = manager->
GetEnergy(p, range, couple);
261 G4cout <<
"G4EmCalculator::GetKinEnergy: Range(mm)= " << range/
mm
262 <<
" KinE(MeV)= " << res/
MeV
282 if(couple && UpdateParticle(p, kinEnergy)) {
283 if(FindEmModel(p, processName, kinEnergy)) {
285 FindLambdaTable(p, processName, kinEnergy);
287 G4VEmProcess* emproc = FindDiscreteProcess(p, processName);
290 }
else if(currentLambda) {
292 res = (((*currentLambda)[idx])->Value(e))*chargeSquare;
298 G4cout <<
"G4EmCalculator::GetXSPerVolume: E(MeV)= " << kinEnergy/
MeV
299 <<
" cross(cm-1)= " << res*
cm
303 G4cout <<
" idx= " << idx <<
" Escaled((MeV)= "
304 << kinEnergy*massRatio
305 <<
" q2= " << chargeSquare;
340 if(x > 0.0) { res = 1.0/
x; }
342 G4cout <<
"G4EmCalculator::GetMeanFreePath: E(MeV)= " << kinEnergy/
MeV
343 <<
" MFP(mm)= " << res/
mm
374 G4cout <<
"### G4EmCalculator: Inverse Range Table for "
391 <<
" in " << currentMaterialName
392 <<
" e(MeV)= " << kinEnergy/
MeV <<
" cut(MeV)= " << cut/
MeV
395 if(UpdateParticle(p, kinEnergy)) {
396 if(FindEmModel(p, processName, kinEnergy)) {
399 if(currentModel->
GetName() ==
"ParamICRU73") {
402 G4cout <<
" ICRU73 ion E(MeV)= " << kinEnergy <<
" ";
404 <<
" DEDX(MeV*cm^2/g)= "
410 G4double escaled = kinEnergy*massRatio;
413 mat, baseParticle, escaled, cut) * chargeSquare;
416 <<
" Escaled(MeV)= " << escaled;
420 if(verbose > 1) {
G4cout <<
" no basePart E(MeV)= " << kinEnergy <<
" "; }
424 <<
" DEDX(MeV*cm^2/g)= "
445 G4cout <<
"At boundary energy(MeV)= " << eth/
MeV
446 <<
" DEDX(MeV/mm)= " << res1*
mm/
MeV
454 if(res1 > 0.0 && escaled > 0.0) {
455 res *= (1.0 + (res0/res1 - 1.0)*eth/escaled);
472 G4cout <<
"After Corrections: DEDX(MeV/mm)= " << res*
mm/
MeV
480 G4cout <<
"Sum: E(MeV)= " << kinEnergy/
MeV
481 <<
" DEDX(MeV/mm)= " << res*
mm/
MeV
483 <<
" cut(MeV)= " << cut/
MeV
485 <<
" in " << currentMaterialName
486 <<
" Zi^2= " << chargeSquare
487 <<
" isIon=" << isIon
503 if(UpdateParticle(part, kinEnergy)) {
506 const std::vector<G4VEnergyLossProcess*> vel =
513 for(
G4int i=0; i<
n; ++i) {
516 if(ActiveForParticle(part, p)) {
519 dedx +=
ComputeDEDX(kinEnergy,part,(vel[i])->GetProcessName(),mat,cut);
536 if(UpdateParticle(part, kinEnergy)) {
539 const std::vector<G4VEnergyLossProcess*> vel =
543 if(mat != cutMaterial) {
553 for(
G4int i=0; i<
n; ++i) {
556 if(ActiveForParticle(part, p)) {
564 dedx +=
ComputeDEDX(kinEnergy,part,(vel[i])->GetProcessName(),
592 G4VEmProcess* nucst = FindDiscreteProcess(p,
"nuclearStopping");
603 <<
" NuclearDEDX(MeV/mm)= " << res*
mm/
MeV
604 <<
" NuclearDEDX(MeV*cm^2/g)= "
622 if(UpdateParticle(p, kinEnergy)) {
623 if(FindEmModel(p, processName, kinEnergy)) {
627 e *= kinEnergy*massRatio;
629 mat, baseParticle, e, aCut, e) * chargeSquare;
634 G4cout <<
"G4EmCalculator::ComputeXSPerVolume: E(MeV)= " << kinEnergy/
MeV
635 <<
" cross(cm-1)= " << res*
cm
636 <<
" cut(keV)= " << aCut/
keV
656 if(UpdateParticle(p, kinEnergy)) {
659 if(FindEmModel(p, processName, kinEnergy)) {
663 e *= kinEnergy*massRatio;
666 baseParticle, e, Z, A, aCut) * chargeSquare;
673 <<
" cross(barn)= " << res/
barn
675 <<
" Z= " << Z <<
" A= " << A/(
g/
mole) <<
" g/mole"
676 <<
" cut(keV)= " << aCut/
keV
693 if(UpdateParticle(p, kinEnergy)) {
695 if(FindEmModel(p, processName, kinEnergy)) {
699 e *= kinEnergy*massRatio;
702 e, aCut) * chargeSquare;
709 <<
" cross(barn)= " << res/
barn
711 <<
" Z= " << Z <<
" shellIdx= " << shellIdx
712 <<
" cut(keV)= " << aCut/
keV
732 if(res > 0.0) { res = 1.0/res; }
765 if(x > 0.0) { mfp = 1.0/
x; }
768 <<
" MFP(mm)= " << mfp/
mm
784 ConvertRangeToEnergy(part, mat, range);
792 if(p != currentParticle) {
796 dynParticle.
SetDefinition(const_cast<G4ParticleDefinition*>(p));
803 currentProcess = FindEnergyLossProcess(p);
804 currentProcessName =
"";
820 && currentParticleName !=
"deuteron"
821 && currentParticleName !=
"triton"
822 && currentParticleName !=
"alpha+"
823 && currentParticleName !=
"helium"
824 && currentParticleName !=
"hydrogen"
828 baseParticle = theGenericIon;
830 G4cout <<
"\n G4EmCalculator::UpdateParticle: isIon 1 "
832 <<
" in " << currentMaterial->
GetName()
833 <<
" e= " << kinEnergy <<
G4endl;
847 G4cout <<
"\n NewIon: massR= "<< massRatio <<
" q2= "
848 << chargeSquare <<
" " << currentProcess <<
G4endl;
860 if(name != currentParticleName) {
863 G4cout <<
"### WARNING: G4EmCalculator::FindParticle fails to find "
884 if(name != currentMaterialName) {
886 if(!currentMaterial) {
887 G4cout <<
"### WARNING: G4EmCalculator::FindMaterial fails to find "
891 return currentMaterial;
899 if(reg !=
"" && reg !=
"world") {
915 if(currentMaterial) {
925 size_t nr = store->size();
927 for(
size_t i=0; i<nr; ++i) {
929 material, ((*store)[i])->GetProductionCuts());
930 if(couple) {
break; }
937 ed <<
"G4EmCalculator::FindCouple: fail for material <"
938 << currentMaterialName <<
">";
939 if(region) { ed <<
" and region " << region->
GetName(); }
940 G4Exception(
"G4EmCalculator::FindCouple",
"em0078",
951 if(!currentMaterial) {
return false; }
952 for (
G4int i=0; i<nLocalMaterials; ++i) {
953 if(material == localMaterials[i] && cut == localCuts[i]) {
954 currentCouple = localCouples[i];
955 currentCoupleIndex = currentCouple->
GetIndex();
961 localMaterials.push_back(material);
962 localCouples.push_back(cc);
963 localCuts.push_back(cut);
966 currentCoupleIndex = currentCouple->
GetIndex();
978 if (!currentLambda || p != lambdaParticle || processName != lambdaName) {
979 lambdaName = processName;
980 currentLambda =
nullptr;
984 if(isIon) { part = theGenericIon; }
987 currentName = processName;
988 currentModel =
nullptr;
997 G4cout <<
"G4VEnergyLossProcess is found out: " << currentName
1001 curProcess = elproc;
1006 G4VEmProcess* proc = FindDiscreteProcess(part, processName);
1010 isApplicable =
true;
1012 G4cout <<
"G4VEmProcess is found out: " << currentName <<
G4endl;
1026 isApplicable =
true;
1028 G4cout <<
"G4VMultipleScattering is found out: " << currentName
1044 isApplicable =
false;
1045 if(!p || !currentMaterial) {
1046 G4cout <<
"G4EmCalculator::FindEmModel WARNING: no particle"
1047 <<
" or materail defined; particle: " << p <<
G4endl;
1048 return isApplicable;
1052 G4double scaledEnergy = kinEnergy*massRatio;
1053 if(isIon) { part = theGenericIon; }
1056 G4cout <<
"## G4EmCalculator::FindEmModel for " << partname
1058 <<
") and " << processName <<
" at E(MeV)= " << scaledEnergy
1060 if(p != part) {
G4cout <<
" GenericIon is the base particle" <<
G4endl; }
1064 currentName = processName;
1065 currentModel =
nullptr;
1066 loweModel =
nullptr;
1077 if(loweModel == currentModel) { loweModel =
nullptr; }
1087 G4VEmProcess* proc = FindDiscreteProcess(part, processName);
1095 if(loweModel == currentModel) { loweModel =
nullptr; }
1109 loweModel =
nullptr;
1113 if(loweModel == currentModel) { loweModel =
nullptr; }
1114 isApplicable =
true;
1128 G4cout <<
" LowEnergy model <" << loweModel->
GetName() <<
">";
1133 return isApplicable;
1146 && currentParticleName !=
"deuteron"
1147 && currentParticleName !=
"triton"
1148 && currentParticleName !=
"He3"
1149 && currentParticleName !=
"alpha"
1150 && currentParticleName !=
"alpha+"
1151 && currentParticleName !=
"helium"
1152 && currentParticleName !=
"hydrogen"
1153 ) { part = theGenericIon; }
1171 const std::vector<G4VEnergyLossProcess*>
v =
1174 for(
G4int i=0; i<
n; ++i) {
1175 if((v[i])->GetProcessName() == processName) {
1177 if(ActiveForParticle(part, p)) {
1193 const std::vector<G4VEmProcess*> v =
1196 for(
G4int i=0; i<
n; ++i) {
1197 if((v[i])->GetProcessName() == processName) {
1199 if(ActiveForParticle(part, p)) {
1215 const std::vector<G4VMultipleScattering*> v =
1218 for(
G4int i=0; i<
n; ++i) {
1219 if((v[i])->GetProcessName() == processName) {
1221 if(ActiveForParticle(part, p)) {
1239 for(
G4int i=0; i<nproc; ++i) {
1240 if(processName == (*pv)[i]->GetProcessName()) {
1258 for(
G4int i=0; i<
n; ++i) {
1259 if((*pv)[i] == proc) {
1272 currentMaterial = mat;
1273 currentMaterialName = mat->
GetName();
1275 currentMaterial = 0;
1276 currentMaterialName =
"";
1289 void G4EmCalculator::CheckMaterial(
G4int Z)
1292 if(currentMaterial) {
1294 for(
size_t i=0; i<
nn; ++i) {
G4PhysicsTable * LambdaTable() const
G4ProductionCuts * GetProductionCuts() const
G4double LowEnergyLimit() const
virtual G4double CrossSectionPerVolume(const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
G4Material * FindOrBuildMaterial(const G4String &name, G4bool isotopes=true, G4bool warning=false)
G4VEmModel * SelectModel(G4double kinEnergy, size_t idx)
G4Region * GetRegion(const G4String &name, G4bool verbose=true) const
virtual void InitialiseForElement(const G4ParticleDefinition *, G4int Z)
const G4String & GetName() const
G4ParticleDefinition * FindParticle(G4int PDGEncoding)
static G4LossTableManager * Instance()
const G4MaterialCutsCouple * FindCouple(const G4Material *, const G4Region *r=nullptr)
static constexpr double mm
virtual G4double GetShellIonisationCrossSectionPerAtom(const G4ParticleDefinition *, G4int Z, G4AtomicShellEnumerator shell, G4double kinE, const G4Material *mat=nullptr)=0
virtual void CorrectionsAlongStep(const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double &eloss, G4double &niel, G4double length)
void SetDynamicMassCharge(G4double massratio, G4double charge2ratio)
std::ostringstream G4ExceptionDescription
G4double GetDEDX(const G4ParticleDefinition *aParticle, G4double kineticEnergy, const G4MaterialCutsCouple *couple)
const G4Material * FindMaterial(const G4String &)
static constexpr double cm2
G4double GetRange(G4double kinEnergy, const G4ParticleDefinition *, const G4Material *, const G4Region *r=nullptr)
G4Material * FindOrBuildSimpleMaterial(G4int Z, G4bool warning=false)
static G4Material * GetMaterial(const G4String &name, G4bool warning=true)
G4double LowestElectronEnergy() const
G4PhysicsTable * RangeTableForLoss() const
G4double EffectiveChargeSquareRatio(const G4ParticleDefinition *, const G4Material *, G4double kineticEnergy)
const G4String & GetName() const
void PrintRangeTable(const G4ParticleDefinition *)
G4ParticleDefinition * GetIon(G4int Z, G4int A, G4int lvl=0)
virtual void SetupForMaterial(const G4ParticleDefinition *, const G4Material *, G4double kineticEnergy)
G4double GetDensity() const
G4double ComputeTotalDEDX(G4double kinEnergy, const G4ParticleDefinition *, const G4Material *, G4double cut=DBL_MAX)
G4double ComputeElectronicDEDX(G4double kinEnergy, const G4ParticleDefinition *, const G4Material *mat, G4double cut=DBL_MAX)
G4double GetCSDARange(G4double kinEnergy, const G4ParticleDefinition *, const G4Material *, const G4Region *r=nullptr)
G4PhysicsTable * GetCrossSectionTable()
G4double EffectiveChargeCorrection(const G4ParticleDefinition *, const G4Material *, G4double kineticEnergy)
void SetupMaterial(const G4Material *)
void PrintInverseRangeTable(const G4ParticleDefinition *)
const G4Element * GetElement(G4int iel) const
const std::vector< G4VEmProcess * > & GetEmProcessVector()
G4double GetKinEnergy(G4double range, const G4ParticleDefinition *, const G4Material *, const G4Region *r=nullptr)
static G4NistManager * Instance()
const G4String & GetParticleName() const
G4double GetCSDARange(const G4ParticleDefinition *aParticle, G4double kineticEnergy, const G4MaterialCutsCouple *couple)
static const G4double reg
static constexpr double gram
static G4RegionStore * GetInstance()
G4bool BuildCSDARange() const
G4double GetShellIonisationCrossSectionPerAtom(const G4String &part, G4int Z, G4AtomicShellEnumerator shell, G4double kinEnergy)
G4PhysicsTable * LambdaTable() const
G4double ComputeShellIonisationCrossSectionPerAtom(const G4String &part, G4int Z, G4AtomicShellEnumerator shell, G4double kinEnergy, const G4Material *mat=nullptr)
G4IonTable * GetIonTable() const
virtual void InitialiseForMaterial(const G4ParticleDefinition *, const G4Material *)
G4GLOB_DLL std::ostream G4cout
double A(double temperature)
G4PhysicsTable * DEDXTable() const
G4double ComputeDEDX(G4double kinEnergy, const G4ParticleDefinition *, const G4String &processName, const G4Material *, G4double cut=DBL_MAX)
const G4ParticleDefinition const G4Material *G4double range
virtual G4double ComputeShellIonisationCrossSectionPerAtom(const G4ParticleDefinition *, G4int Z, G4AtomicShellEnumerator shell, G4double kinE, const G4Material *mat=nullptr)=0
G4EmCorrections * EmCorrections()
G4double GetRangeFromRestricteDEDX(G4double kinEnergy, const G4ParticleDefinition *, const G4Material *, const G4Region *r=nullptr)
G4double GetDEDX(G4double kinEnergy, const G4ParticleDefinition *, const G4Material *, const G4Region *r=nullptr)
G4double ComputeNuclearDEDX(G4double kinEnergy, const G4ParticleDefinition *, const G4Material *)
static constexpr double cm
const G4String & GetParticleType() const
const G4ParticleDefinition * BaseParticle() const
G4double GetEnergy(const G4ParticleDefinition *aParticle, G4double range, const G4MaterialCutsCouple *couple)
const G4String & GetProcessName() const
const std::vector< G4VEnergyLossProcess * > & GetEnergyLossProcessVector()
virtual G4double ComputeCrossSectionPerAtom(const G4ParticleDefinition *, G4double kinEnergy, G4double Z, G4double A=0., G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
void SetKineticEnergy(G4double aEnergy)
G4double GetRangeFromRestricteDEDX(const G4ParticleDefinition *aParticle, G4double kineticEnergy, const G4MaterialCutsCouple *couple)
static constexpr double eV
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
G4double GetCrossSectionPerVolume(G4double kinEnergy, const G4ParticleDefinition *, const G4String &processName, const G4Material *, const G4Region *r=nullptr)
virtual G4double GetChargeSquareRatio(const G4ParticleDefinition *, const G4Material *, G4double kineticEnergy)
G4VEnergyLossProcess * GetEnergyLossProcess(const G4ParticleDefinition *)
static G4ProductionCutsTable * GetProductionCutsTable()
void SetFluctuationFlag(G4bool val)
static G4Positron * Positron()
static G4GenericIon * GenericIon()
virtual G4double ComputeDEDXPerVolume(const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=DBL_MAX)
G4VEmModel * SelectModelForMaterial(G4double kinEnergy, size_t &idxRegion) const
const G4ParticleDefinition * FindParticle(const G4String &)
G4PhysicsTable * InverseRangeTable() const
static constexpr double nm
G4double GetPDGMass() const
const G4MaterialCutsCouple * GetMaterialCutsCouple(G4int i) const
static G4ParticleTable * GetParticleTable()
T max(const T t1, const T t2)
brief Return the largest of the two arguments
G4double ComputeMeanFreePath(G4double kinEnergy, const G4ParticleDefinition *, const G4String &processName, const G4Material *, G4double cut=0.0)
G4ProcessManager * GetProcessManager() const
G4double ComputeCrossSectionPerShell(G4double kinEnergy, const G4ParticleDefinition *, const G4String &processName, G4int Z, G4int shellIdx, G4double cut=0.0)
const G4Region * FindRegion(const G4String &)
const G4ParticleDefinition * FindIon(G4int Z, G4int A)
static G4EmParameters * Instance()
G4double ComputeCrossSectionPerVolume(G4double kinEnergy, const G4ParticleDefinition *, const G4String &processName, const G4Material *, G4double cut=0.0)
static constexpr double GeV
virtual G4double ComputeCrossSectionPerShell(const G4ParticleDefinition *, G4int Z, G4int shellIdx, G4double kinEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
G4double ComputeGammaAttenuationLength(G4double kinEnergy, const G4Material *)
G4bool GetProcessActivation(G4VProcess *aProcess) const
static G4Electron * Electron()
static constexpr double MeV
G4double ComputeCrossSectionPerAtom(G4double kinEnergy, const G4ParticleDefinition *, const G4String &processName, G4double Z, G4double A, G4double cut=0.0)
const G4String & GetName() const
G4VProcess * FindProcess(const G4ParticleDefinition *part, const G4String &processName)
size_t GetNumberOfElements() const
void SetVerbose(G4int val)
G4VEmModel * GetModelByIndex(G4int idx=0, G4bool ver=false) const
void PrintDEDXTable(const G4ParticleDefinition *)
G4VAtomDeexcitation * AtomDeexcitation()
void SetDefinition(const G4ParticleDefinition *aParticleDefinition)
G4double ComputeEnergyCutFromRangeCut(G4double range, const G4ParticleDefinition *, const G4Material *)
static constexpr double barn
G4double GetPDGCharge() const
static constexpr double keV
G4double CrossSectionPerVolume(G4double kineticEnergy, const G4MaterialCutsCouple *couple)
static constexpr double mole
G4double ComputeDEDXForCutInRange(G4double kinEnergy, const G4ParticleDefinition *, const G4Material *mat, G4double rangecut=DBL_MAX)
G4double GetMeanFreePath(G4double kinEnergy, const G4ParticleDefinition *, const G4String &processName, const G4Material *, const G4Region *r=nullptr)
const std::vector< G4VMultipleScattering * > & GetMultipleScatteringVector()
G4VEmModel * SelectModelForMaterial(G4double kinEnergy, size_t &idx) const
G4ProcessVector * GetProcessList() const