lars/geant4
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Import Geant4 11.1.1 source tree

Browse Source
This commit is contained in:
Gabriele Cosmo
2023-02-14 13:57:32 +01:00
parent 9f34590941
commit 84a556a9dc
312 changed files with 35018 additions and 36005 deletions
Download Patch File Download Diff File Expand all files Collapse all files
source/processes/hadronic/processes/History
+18
View File
@@ -6,6 +6,24 @@ It must **not** be used as a substitute for writing good git commit messages!
-------------------------------------------------------------------------------
## 2023-02-10 Vladimir Ivanchenko (hadr-proc-V11-00-13)
- G4NeutronGeneralProcess - added extra Set method
## 2023-02-01 Vladimir Ivanchenko (hadr-proc-V11-00-12)
- G4HadronElasticProcess, G4NeutronGeneralProcess - fixed problems identified
by Coverity
## 2023-01-17 Vladimir Ivanchenko
- G4HadronElasticProcess - removed obsolete (try/catch construction not
applicable for hadron elastic models; added forgotten integral approach
for charged particles
## 2022-12-30 Vladimir Ivanchenko
- G4NeutronGeneralProcess - optimized code: initialisation methods are
moved to the source, avoid double instantiation of capture cross section;
reduced number of calls to cross section, added cross section data store
pointer
## 2022-11-26 Gabriele Cosmo (hadr-proc-V11-00-11)
- Fixed compilation warnings for implicit type conversions on macOS/XCode 14.1.
source/processes/hadronic/processes/include/G4NeutronGeneralProcess.hh
+29 -47
View File
@@ -57,6 +57,7 @@ class G4Track;
class G4ParticleDefinition;
class G4VParticleChange;
class G4VCrossSectionDataSet;
class G4CrossSectionDataStore;
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
@@ -102,14 +103,16 @@ public:
// Temporary method
G4int GetSubProcessSubType() const;
void SetInelasticProcess(G4HadronicProcess*);
void SetElasticProcess(G4HadronicProcess*);
void SetCaptureProcess(G4HadronicProcess*);
inline const G4VProcess* GetSelectedProcess() const;
inline void SetInelasticProcess(G4HadronicProcess*);
inline void SetElasticProcess(G4HadronicProcess*);
inline void SetCaptureProcess(G4HadronicProcess*);
inline void SetTimeLimit(G4double val);
inline void SetMinEnergyLimit(G4double val);
// hide copy constructor and assignment operator
G4NeutronGeneralProcess(G4NeutronGeneralProcess &) = delete;
G4NeutronGeneralProcess & operator=
@@ -125,15 +128,15 @@ protected:
inline G4double GetProbability(size_t idxt);
inline void SelectedProcess(const G4Step& step, G4HadronicProcess* ptr,
G4VCrossSectionDataSet* xs);
void SelectHadProcess(const G4Track&, const G4Step&, G4HadronicProcess*);
G4CrossSectionDataStore*);
private:
// partial cross section
G4double ComputeCrossSection(G4VCrossSectionDataSet*, const G4Material*,
G4double kinEnergy, G4double loge);
G4double kinEnergy, G4double loge);
G4VCrossSectionDataSet* InitialisationXS(G4HadronicProcess*);
// total cross section
inline void CurrentCrossSection(const G4Track&);
@@ -147,10 +150,14 @@ private:
G4HadronicProcess* fCapture = nullptr;
G4HadronicProcess* fSelectedProc = nullptr;
G4VCrossSectionDataSet* fInelasticXS;
G4VCrossSectionDataSet* fElasticXS;
G4VCrossSectionDataSet* fCaptureXS;
G4VCrossSectionDataSet* fXS = nullptr;
G4VCrossSectionDataSet* fInelasticXS = nullptr;
G4VCrossSectionDataSet* fElasticXS = nullptr;
G4VCrossSectionDataSet* fCaptureXS = nullptr;
G4CrossSectionDataStore* fXSSInelastic = nullptr;
G4CrossSectionDataStore* fXSSElastic = nullptr;
G4CrossSectionDataStore* fXSSCapture = nullptr;
G4CrossSectionDataStore* fCurrentXSS = nullptr;
const G4ParticleDefinition* fNeutron;
const G4Material* fCurrMat = nullptr;
@@ -178,31 +185,6 @@ private:
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
inline void
G4NeutronGeneralProcess::SetInelasticProcess(G4HadronicProcess* ptr)
{
fInelastic = ptr;
ptr->AddDataSet(fInelasticXS);
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
inline void G4NeutronGeneralProcess::SetElasticProcess(G4HadronicProcess* ptr)
{
fElastic = ptr;
ptr->AddDataSet(fElasticXS);
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
inline void G4NeutronGeneralProcess::SetCaptureProcess(G4HadronicProcess* ptr)
{
fCapture = ptr;
ptr->AddDataSet(fCaptureXS);
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
inline G4double
G4NeutronGeneralProcess::ComputeGeneralLambda(std::size_t idxe, std::size_t idxt)
{
@@ -224,11 +206,11 @@ inline G4double G4NeutronGeneralProcess::GetProbability(std::size_t idxt)
inline void
G4NeutronGeneralProcess::SelectedProcess(const G4Step& step,
G4HadronicProcess* ptr,
G4VCrossSectionDataSet* xs)
G4CrossSectionDataStore* xs)
{
fSelectedProc = ptr;
fXS = xs;
fCurrentXSS = xs;
step.GetPostStepPoint()->SetProcessDefinedStep(ptr);
}
@@ -245,18 +227,11 @@ inline void G4NeutronGeneralProcess::CurrentCrossSection(const G4Track& track)
{
G4double energy = track.GetKineticEnergy();
const G4Material* mat = track.GetMaterial();
G4bool recompute = false;
if(mat != fCurrMat) {
if(mat != fCurrMat || energy != fCurrE) {
fCurrMat = mat;
matIndex = mat->GetIndex();
recompute = true;
}
if(energy != fCurrE) {
fCurrE = energy;
fCurrLogE = track.GetDynamicParticle()->GetLogKineticEnergy();
recompute = true;
}
if(recompute) {
fLambda = (energy <= fMiddleEnergy) ? ComputeGeneralLambda(0, 0)
: ComputeGeneralLambda(1, 3);
currentInteractionLength = 1.0/fLambda;
@@ -272,4 +247,11 @@ inline void G4NeutronGeneralProcess::SetTimeLimit(G4double val)
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
inline void G4NeutronGeneralProcess::SetMinEnergyLimit(G4double val)
{
fMinEnergy = val;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
#endif
source/processes/hadronic/processes/src/G4HadronElasticProcess.cc
+38 -47
View File
@@ -68,30 +68,41 @@ G4HadronElasticProcess::PostStepDoIt(const G4Track& track,
theTotalResult->ProposeWeight(weight);
// For elastic scattering, _any_ result is considered an interaction
ClearNumberOfInteractionLengthLeft();
theNumberOfInteractionLengthLeft = -1.0;
G4double kineticEnergy = track.GetKineticEnergy();
const G4DynamicParticle* dynParticle = track.GetDynamicParticle();
G4double kineticEnergy = dynParticle->GetKineticEnergy();
G4TrackStatus status = track.GetTrackStatus();
if(kineticEnergy == 0.0 || track.GetTrackStatus() != fAlive) {
return theTotalResult;
}
const G4DynamicParticle* dynParticle = track.GetDynamicParticle();
const G4ParticleDefinition* part = dynParticle->GetDefinition();
const G4Material* material = track.GetMaterial();
// check only for charged particles
if(fXSType != fHadNoIntegral) {
mfpKinEnergy = DBL_MAX;
G4double xs = aScaleFactor*
theCrossSectionDataStore->ComputeCrossSection(dynParticle, material);
if(xs < theLastCrossSection*G4UniformRand()) {
// No interaction
return theTotalResult;
}
}
const G4ParticleDefinition* part = dynParticle->GetDefinition();
G4Nucleus* targNucleus = GetTargetNucleusPointer();
// Select element
const G4Element* elm =
GetCrossSectionDataStore()->SampleZandA(dynParticle, material, *targNucleus);
theCrossSectionDataStore->SampleZandA(dynParticle, material, *targNucleus);
// Initialize the hadronic projectile from the track
G4HadProjectile theProj(track);
G4HadronicInteraction* hadi = nullptr;
G4HadFinalState* result = nullptr;
if(fDiffraction)
{
if(nullptr != fDiffraction) {
G4double ratio =
fDiffractionRatio->ComputeRatio(part, kineticEnergy,
targNucleus->GetZ_asInt(),
@@ -108,7 +119,8 @@ G4HadronElasticProcess::PostStepDoIt(const G4Track& track,
G4ExceptionDescription ed;
aR.Report(ed);
ed << "Call for " << fDiffraction->GetModelName() << G4endl;
ed << "Target element "<< elm->GetName()<<" Z= "
ed << part->GetParticleName()
<< " off target element " << elm->GetName() << " Z= "
<< targNucleus->GetZ_asInt()
<< " A= " << targNucleus->GetA_asInt() << G4endl;
DumpState(track,"ApplyYourself",ed);
@@ -118,9 +130,7 @@ G4HadronElasticProcess::PostStepDoIt(const G4Track& track,
}
// Check the result for catastrophic energy non-conservation
result = CheckResult(theProj, *targNucleus, result);
result->SetTrafoToLab(theProj.GetTrafoToLab());
ClearNumberOfInteractionLengthLeft();
// The following method of the base class takes care also of setting
// the creator model ID for the secondaries that are created
@@ -134,28 +144,25 @@ G4HadronElasticProcess::PostStepDoIt(const G4Track& track,
}
// ordinary elastic scattering
try
{
hadi = ChooseHadronicInteraction( theProj, *targNucleus, material, elm );
}
catch(G4HadronicException & aE)
{
G4ExceptionDescription ed;
aE.Report(ed);
ed << "Target element "<< elm->GetName()<<" Z= "
<< targNucleus->GetZ_asInt() << " A= "
<< targNucleus->GetA_asInt() << G4endl;
DumpState(track,"ChooseHadronicInteraction",ed);
ed << " No HadronicInteraction found out" << G4endl;
G4Exception("G4HadronElasticProcess::PostStepDoIt", "had005",
FatalException, ed);
}
hadi = ChooseHadronicInteraction( theProj, *targNucleus, material, elm );
if(nullptr == hadi) {
G4ExceptionDescription ed;
ed << part->GetParticleName()
<< " off target element " << elm->GetName() << " Z= "
<< targNucleus->GetZ_asInt() << " A= "
<< targNucleus->GetA_asInt() << G4endl;
DumpState(track,"ChooseHadronicInteraction",ed);
ed << " No HadronicInteraction found out" << G4endl;
G4Exception("G4HadronElasticProcess::PostStepDoIt", "had005",
FatalException, ed);
return theTotalResult;
}
size_t idx = track.GetMaterialCutsCouple()->GetIndex();
G4double tcut = (*(G4ProductionCutsTable::GetProductionCutsTable()
->GetEnergyCutsVector(3)))[idx];
hadi->SetRecoilEnergyThreshold(tcut);
/*
if(verboseLevel>1) {
G4cout << "G4HadronElasticProcess::PostStepDoIt for "
<< part->GetParticleName()
@@ -164,24 +171,8 @@ G4HadronElasticProcess::PostStepDoIt(const G4Track& track,
<< " A= " << targNucleus->GetA_asInt()
<< " Tcut(MeV)= " << tcut << G4endl;
}
try
{
result = hadi->ApplyYourself( theProj, *targNucleus);
}
catch(G4HadronicException & aR)
{
G4ExceptionDescription ed;
aR.Report(ed);
ed << "Call for " << hadi->GetModelName() << G4endl;
ed << "Target element "<< elm->GetName()<<" Z= "
<< targNucleus->GetZ_asInt()
<< " A= " << targNucleus->GetA_asInt() << G4endl;
DumpState(track,"ApplyYourself",ed);
ed << " ApplyYourself failed" << G4endl;
G4Exception("G4HadronElasticProcess::PostStepDoIt", "had006",
FatalException, ed);
}
*/
result = hadi->ApplyYourself( theProj, *targNucleus);
// Check the result for catastrophic energy non-conservation
// cannot be applied because is not guranteed that recoil
@@ -191,7 +182,7 @@ G4HadronElasticProcess::PostStepDoIt(const G4Track& track,
// directions
G4ThreeVector indir = track.GetMomentumDirection();
G4ThreeVector outdir = result->GetMomentumChange();
/*
if(verboseLevel>1) {
G4cout << "Efin= " << result->GetEnergyChange()
<< " de= " << result->GetLocalEnergyDeposit()
@@ -199,7 +190,7 @@ G4HadronElasticProcess::PostStepDoIt(const G4Track& track,
<< " dir= " << outdir
<< G4endl;
}
*/
// energies
G4double edep = std::max(result->GetLocalEnergyDeposit(), 0.0);
G4double efinal = std::max(result->GetEnergyChange(), 0.0);
source/processes/hadronic/processes/src/G4NeutronGeneralProcess.cc
+66 -43
View File
@@ -48,8 +48,8 @@
#include "G4NeutronGeneralProcess.hh"
#include "G4PhysicalConstants.hh"
#include "G4SystemOfUnits.hh"
#include "G4ProcessManager.hh"
#include "G4HadronicProcess.hh"
#include "G4CrossSectionDataStore.hh"
#include "G4Step.hh"
#include "G4Track.hh"
#include "G4ParticleDefinition.hh"
@@ -62,7 +62,6 @@
#include "G4MaterialTable.hh"
#include "G4Element.hh"
#include "G4Neutron.hh"
#include "G4Nucleus.hh"
#include "G4NeutronInelasticXS.hh"
#include "G4NeutronElasticXS.hh"
#include "G4NeutronCaptureXS.hh"
@@ -87,14 +86,6 @@ G4NeutronGeneralProcess::G4NeutronGeneralProcess(const G4String& pname)
SetVerboseLevel(1);
SetProcessSubType(fNeutronGeneral);
fElasticXS = new G4NeutronElasticXS();
fInelasticXS = new G4NeutronInelasticXS();
fCaptureXS = new G4NeutronCaptureXS();
AddDataSet(fElasticXS);
AddDataSet(fInelasticXS);
AddDataSet(fCaptureXS);
fNeutron = G4Neutron::Neutron();
if(G4Threading::IsWorkerThread()) {
@@ -121,6 +112,58 @@ G4bool G4NeutronGeneralProcess::IsApplicable(const G4ParticleDefinition&)
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
void G4NeutronGeneralProcess::SetInelasticProcess(G4HadronicProcess* ptr)
{
fInelastic = ptr;
fXSSInelastic = ptr->GetCrossSectionDataStore();
fInelasticXS = InitialisationXS(ptr);
if(nullptr == fInelasticXS) {
fInelasticXS = new G4NeutronInelasticXS();
ptr->AddDataSet(fInelasticXS);
}
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
void G4NeutronGeneralProcess::SetElasticProcess(G4HadronicProcess* ptr)
{
fElastic = ptr;
fXSSElastic = ptr->GetCrossSectionDataStore();
fElasticXS = InitialisationXS(ptr);
if(nullptr == fElasticXS) {
fElasticXS = new G4NeutronElasticXS();
ptr->AddDataSet(fElasticXS);
}
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
void G4NeutronGeneralProcess::SetCaptureProcess(G4HadronicProcess* ptr)
{
fCapture = ptr;
fXSSCapture = ptr->GetCrossSectionDataStore();
fCaptureXS = InitialisationXS(ptr);
if(nullptr == fCaptureXS) {
fCaptureXS = new G4NeutronCaptureXS();
ptr->AddDataSet(fCaptureXS);
}
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
G4VCrossSectionDataSet*
G4NeutronGeneralProcess::InitialisationXS(G4HadronicProcess* proc)
{
G4VCrossSectionDataSet* ptr = nullptr;
auto xsv = proc->GetCrossSectionDataStore()->GetDataSetList();
if(!xsv.empty()) {
ptr = xsv[0];
}
return ptr;
}
//....Ooooo0ooooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
void G4NeutronGeneralProcess::PreparePhysicsTable(const G4ParticleDefinition& part)
{
if(1 < verboseLevel) {
@@ -211,7 +254,6 @@ void G4NeutronGeneralProcess::BuildPhysicsTable(const G4ParticleDefinition& part
fElastic->BuildPhysicsTable(part);
fInelastic->BuildPhysicsTable(part);
fCapture->BuildPhysicsTable(part);
fCaptureXS->BuildPhysicsTable(part);
if(isMaster) {
std::size_t nmat = G4Material::GetNumberOfMaterials();
@@ -349,7 +391,7 @@ G4double G4NeutronGeneralProcess::PostStepGetPhysicalInteractionLength(
G4VParticleChange* G4NeutronGeneralProcess::PostStepDoIt(const G4Track& track,
const G4Step& step)
{
fSelectedProc = nullptr;
fSelectedProc = this;
// time limit
if(0.0 == fLambda) {
theTotalResult->Initialize(track);
@@ -365,46 +407,26 @@ G4VParticleChange* G4NeutronGeneralProcess::PostStepDoIt(const G4Track& track,
*/
if (0 == idxEnergy) {
if(q <= GetProbability(1)) {
SelectedProcess(step, fElastic, fElasticXS);
SelectedProcess(step, fElastic, fXSSElastic);
} else if(q <= GetProbability(2)) {
SelectedProcess(step, fInelastic, fInelasticXS);
SelectedProcess(step, fInelastic, fXSSInelastic);
} else {
SelectedProcess(step, fCapture, fCaptureXS);
SelectedProcess(step, fCapture, fXSSCapture);
}
} else {
if(q <= GetProbability(4)) {
SelectedProcess(step, fInelastic, fInelasticXS);
SelectedProcess(step, fInelastic, fXSSInelastic);
} else {
SelectedProcess(step, fElastic, fElasticXS);
SelectedProcess(step, fElastic, fXSSElastic);
}
}
const G4Element* elm = fCurrMat->GetElement(0);
G4int nelm = (G4int)fCurrMat->GetNumberOfElements();
if(1 < nelm) {
auto natom = fCurrMat->GetVecNbOfAtomsPerVolume();
G4double sig = 0.0;
for(G4int i=0; i<nelm; ++i) {
sig += natom[i] *
fXS->ComputeCrossSectionPerElement(fCurrE, fCurrLogE, fNeutron,
fCurrMat->GetElement(i),
fCurrMat);
fXsec[i] = sig;
}
sig *= G4UniformRand();
for(G4int i=0; i<nelm; ++i) {
if(fXsec[i] >= sig) {
elm = fCurrMat->GetElement(i);
break;
}
}
// total cross section is needed for selection of an element
if(fCurrMat->GetNumberOfElements() > 1) {
fCurrentXSS->ComputeCrossSection(track.GetDynamicParticle(), fCurrMat);
}
fSelectedProc->GetCrossSectionDataStore()->SetForcedElement(elm);
const G4Isotope* iso = fXS->SelectIsotope(elm, fCurrE, fCurrLogE);
fSelectedProc->GetTargetNucleusPointer()->SetIsotope(iso);
/*
G4cout << "## neutron E(MeV)=" << fCurrE << " "
G4cout << "## neutron E(MeV)=" << fCurrE << " inside " << fCurrMat->GetName()
<< fSelectedProc->GetProcessName()
<< " on Z=" << iso->GetZ() << " A=" << iso->GetN()
<< " time(ns)=" << track.GetGlobalTime()/ns << G4endl;
*/
// sample secondaries
@@ -457,7 +479,7 @@ void G4NeutronGeneralProcess::ProcessDescription(std::ostream& out) const
const G4String& G4NeutronGeneralProcess::GetSubProcessName() const
{
return (fSelectedProc) ? fSelectedProc->GetProcessName()
return (nullptr != fSelectedProc) ? fSelectedProc->GetProcessName()
: G4VProcess::GetProcessName();
}
@@ -465,7 +487,8 @@ const G4String& G4NeutronGeneralProcess::GetSubProcessName() const
G4int G4NeutronGeneralProcess::GetSubProcessSubType() const
{
return (fSelectedProc) ? fSelectedProc->GetProcessSubType() : 16;
return (nullptr != fSelectedProc) ? fSelectedProc->GetProcessSubType()
: fNeutronGeneral;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....