Skip to content

Commit

Permalink
Merge pull request #59 from civanch/geant4_composits
Browse files Browse the repository at this point in the history 
Fixed Geant4 compound materials construction
  • Loading branch information
@smuzaffar
smuzaffar committed Apr 15, 2022
2 parents 361837f + 7a2b826 commit 48e04df
Show file tree
Hide file tree
Showing 7 changed files with 110 additions and 130 deletions.
13 changes: 13 additions & 0 deletions source/materials/History
View file
Expand Up @@ -16,6 +16,19 @@ committal in the CVS repository !
* Reverse chronological order (last date on top), please *
----------------------------------------------------------

## 2022-04-14 Vladimir Ivanchenko (materials-V10-07-28)
## 2022-04-12 Vladimir Ivanchenko
- G4DensityEffectData, G4IonisParamMat - fixed density effect correction
for liquid hydrogen (problem #2346), code clean-up
- G4NistMaterialBuilder - fixed mean ionisation potentilal of carbon
to the current NIST value from 81 eV to 78 eV (problem #2474)
- G4Material - reorganized AddElement(..) and AddMaterial(..) methods,
temporary vector are used that allows addition of the same elements
during initialisation, when addition completed the vector of elements
and arrays of fractions are created, filled and temporary vectors
are deleted, as a result the code become more transparent and
fixing problem #2486

04-02-22 V.Ivanchenko (materials-V10-07-27)
- G4SandiaTable - fixed compilation warning on unused variable

Expand Down
2 changes: 1 addition & 1 deletion source/materials/include/G4DensityEffectData.hh
View file
Expand Up @@ -61,7 +61,7 @@ public:
~G4DensityEffectData();

// return index by Z, -1 if material is not in the table
G4int GetElementIndex(G4int Z, G4State mState) const;
G4int GetElementIndex(G4int Z, G4State st = kStateUndefined) const;

// return index by material name, -1 if material is not in the table
G4int GetIndex(const G4String& matName) const;
Expand Down
5 changes: 5 additions & 0 deletions source/materials/include/G4Material.hh
View file
Expand Up @@ -351,6 +351,11 @@ private:
G4int fIdxComponent; // Index of a new component
G4bool fMassFraction; // Flag of the method to add components

// For composites built
std::vector<G4int>* fAtoms;
std::vector<G4double>* fElmFrac;
std::vector<const G4Element*>* fElm;

// For composites built via AddMaterial()
std::map<G4Material*, G4double> fMatComponents;

Expand Down
8 changes: 2 additions & 6 deletions source/materials/src/G4DensityEffectData.cc
View file
Expand Up @@ -1294,13 +1294,9 @@ void G4DensityEffectData::Initialize()
AddMaterial(M277,"G4_GRAPHITE_POROUS");
}

G4int G4DensityEffectData::GetElementIndex(G4int Z, G4State st) const
G4int G4DensityEffectData::GetElementIndex(G4int Z, G4State) const
{
G4int idx = -1;
if(Z > 0 && Z < NDENSELEM) {
if(st == state[Z] || st == kStateUndefined) { idx = indexZ[Z]; }
}
return idx;
return (Z >= 0 && Z < NDENSELEM) ? indexZ[Z] : -1;
}

G4int G4DensityEffectData::GetIndex(const G4String& matName) const
Expand Down
5 changes: 3 additions & 2 deletions source/materials/src/G4IonisParamMat.cc
View file
Expand Up @@ -212,11 +212,12 @@ void G4IonisParamMat::ComputeDensityEffectParameters()
G4double corr = 0.0;

if(idx < 0 && 1 == nelm) {
idx = fDensityData->GetElementIndex(Z0, fMaterial->GetState());
G4int z = (1 == Z0 && State == kStateLiquid) ? 0 : Z0;
idx = fDensityData->GetElementIndex(z);

// Correction for base material or for non-nominal density
// Except cases of very different density defined in user code
if(idx >= 0) {
if(idx >= 0 && 0 < z) {
G4double dens = nist->GetNominalDensity(Z0);
if(dens <= 0.0) { idx = -1; }
else {
Expand Down
205 changes: 85 additions & 120 deletions source/materials/src/G4Material.cc
View file
Expand Up @@ -167,15 +167,6 @@ G4Material::G4Material(const G4String& name, G4double density,
fPressure = pressure;

fNbComponents = nComponents;
theElementVector = new G4ElementVector();
theElementVector->reserve(fNbComponents);

fAtomsVector = new G4int[fNbComponents];
fMassFractionVector = new G4double[fNbComponents];
for(G4int i=0; i<fNbComponents; ++i) {
fAtomsVector[i] = 0;
fMassFractionVector[i] = 0.0;
}
fMassFraction = true;

if (fState == kStateUndefined)
Expand Down Expand Up @@ -301,7 +292,7 @@ void G4Material::ComputeDerivedQuantities()
// Number of atoms per volume (per element), total nb of electrons per volume
G4double Zi, Ai;
fTotNbOfAtomsPerVolume = 0.;
if (fVecNbOfAtomsPerVolume) { delete [] fVecNbOfAtomsPerVolume; }
delete [] fVecNbOfAtomsPerVolume;
fVecNbOfAtomsPerVolume = new G4double[fNumberOfElements];
fTotNbOfElectPerVolume = 0.;
fFreeElecDensity = 0.;
Expand Down Expand Up @@ -369,7 +360,11 @@ void
G4Material::AddElementByNumberOfAtoms(const G4Element* elm, G4int nAtoms)
{
// perform checks consistency
if(0 == fIdxComponent) { fMassFraction = false; }
if(0 == fIdxComponent) {
fMassFraction = false;
fAtoms = new std::vector<G4int>;
fElm = new std::vector<const G4Element*>;
}
if(fIdxComponent >= fNbComponents) {
G4ExceptionDescription ed;
ed << "For material " << fName << " and added element "
Expand All @@ -381,51 +376,51 @@ G4Material::AddElementByNumberOfAtoms(const G4Element* elm, G4int nAtoms)
}
if(fMassFraction) {
G4ExceptionDescription ed;
G4cout << "For material " << fName << " and added element "
<< elm->GetName() << " with Natoms=" << nAtoms
<< " problem: cannot add by number of atoms after "
<< "addition of elements by mass fraction";
ed << "For material " << fName << " and added element "
<< elm->GetName() << " with Natoms=" << nAtoms
<< " problem: cannot add by number of atoms after "
<< "addition of elements by mass fraction";
G4Exception ("G4Material::AddElementByNumberOfAtoms()", "mat031",
FatalException, ed, "");
}
G4Element* element = const_cast<G4Element*>(elm);
// filling
if (fIdxComponent < fNbComponents) {
G4bool isAdded = false;
// filling
G4bool isAdded = false;
if(!fElm->empty()) {
for (G4int i=0; i<fNumberOfElements; ++i) {
if ( element == (*theElementVector)[i] ) {
G4ExceptionDescription ed;
ed << "For material " << fName << " and added element "
<< elm->GetName() << ", Natoms=" << nAtoms
<< ", fIdxComponent=" << fIdxComponent
<< " problem: attempt to add the same element, which already is at idx="
<< i << " with the Natoms=" << fAtomsVector[i];
G4Exception ("G4Material::AddElementByNumberOfAtoms()", "mat031",
JustWarning, ed, "");
fAtomsVector[i] += nAtoms;
if ( elm == (*fElm)[i] ) {
(*fAtoms)[i] += nAtoms;
isAdded = true;
break;
}
}
if(!isAdded) {
theElementVector->push_back(element);
fAtomsVector[fNumberOfElements] = nAtoms;
++fNumberOfElements;
}
}
if(!isAdded) {
fElm->push_back(elm);
fAtoms->push_back(nAtoms);
++fNumberOfElements;
}
++fIdxComponent;
// is filled

// is filled - complete composition of atoms
if (fIdxComponent == fNbComponents) {
// compute proportion by mass
theElementVector = new G4ElementVector();
theElementVector->reserve(fNumberOfElements);
fAtomsVector = new G4int[fNumberOfElements];
fMassFractionVector = new G4double[fNumberOfElements];

G4double Amol = 0.;
for (G4int i=0; i<fNumberOfElements; ++i) {
G4double w = fAtomsVector[i]*(*theElementVector)[i]->GetA();
theElementVector->push_back((*fElm)[i]);
fAtomsVector[i] = (*fAtoms)[i];
G4double w = fAtomsVector[i]*(*fElm)[i]->GetA();
Amol += w;
fMassFractionVector[i] = w;
}
for (G4int i=0; i<fNumberOfElements; ++i) {
fMassFractionVector[i] /= Amol;
}

delete fAtoms;
delete fElm;
fMassOfMolecule = Amol/Avogadro;
ComputeDerivedQuantities();
}
Expand Down Expand Up @@ -462,36 +457,31 @@ G4Material::AddElementByMassFraction(const G4Element* elm, G4double fraction)
ed << "For material " << fName << " and added element "
<< elm->GetName() << ", massFraction= " << fraction
<< ", fIdxComponent=" << fIdxComponent
<< " problem: attempt to add more than the declared number of elements "
<< "; attempt to add more than the declared number of components "
<< fIdxComponent << " >= " << fNbComponents;
G4Exception ("G4Material::AddElementByMassFraction()", "mat031",
FatalException, ed, "");
}
G4Element* element = const_cast<G4Element*>(elm);
if(0 == fIdxComponent) {
fElmFrac = new std::vector<G4double>;
fElm = new std::vector<const G4Element*>;
}

// filling
if (fIdxComponent < fNbComponents) {
G4bool isAdded = false;
G4bool isAdded = false;
if(!fElm->empty()) {
for (G4int i=0; i<fNumberOfElements; ++i) {
if ( element == (*theElementVector)[i] ) {
G4ExceptionDescription ed;
ed << "For material " << fName << " and added element "
<< elm->GetName() << ", massFraction= " << fraction
<< ", fIdxComponent=" << fIdxComponent
<< " problem: attempt to add the same element, which is already at idx="
<< i << " with the fraction " << fMassFractionVector[i];
G4Exception ("G4Material::AddElementByMassFraction()", "mat031",
JustWarning, ed, "");
fMassFractionVector[i] += fraction;
if ( elm == (*fElm)[i] ) {
(*fElmFrac)[i] += fraction;
isAdded = true;
break;
}
}
if(!isAdded) {
theElementVector->push_back(element);
fMassFractionVector[fNumberOfElements] = fraction;
++fNumberOfElements;
}
}
if(!isAdded) {
fElm->push_back(elm);
fElmFrac->push_back(fraction);
++fNumberOfElements;
}
++fIdxComponent;

Expand Down Expand Up @@ -527,15 +517,40 @@ void G4Material::AddMaterial(G4Material* material, G4double fraction)
G4ExceptionDescription ed;
ed << "For material " << fName << " and added material "
<< material->GetName() << ", massFraction= " << fraction
<< " attempt to add more than the declared number of elements "
<< "; attempt to add more than the declared number of components "
<< fIdxComponent << " >= " << fNbComponents;
G4Exception ("G4Material::AddMaterial()", "mat031", FatalException,
ed, "");
}

if(0 == fIdxComponent) {
fElmFrac = new std::vector<G4double>;
fElm = new std::vector<const G4Element*>;
}

// filling
if (fIdxComponent < fNbComponents) {
fMatComponents[material] = fraction;
G4int nelm = material->GetNumberOfElements();
for(G4int j=0; j<nelm; ++j) {
auto elm = material->GetElement(j);
auto frac = material->GetFractionVector();
G4bool isAdded = false;
if(!fElm->empty()) {
for (G4int i=0; i<fNumberOfElements; ++i) {
if ( elm == (*fElm)[i] ) {
(*fElmFrac)[i] += fraction*frac[j];
isAdded = true;
break;
}
}
}
if(!isAdded) {
fElm->push_back(elm);
fElmFrac->push_back(fraction*frac[j]);
++fNumberOfElements;
}
}

fMatComponents[material] = fraction;
++fIdxComponent;

// is filled
Expand All @@ -547,70 +562,20 @@ void G4Material::AddMaterial(G4Material* material, G4double fraction)
void G4Material::FillVectors()
{
// there are material components
if(!fMatComponents.empty()) {
G4int nel = fNumberOfElements;
// check list of materials
for(auto & x : fMatComponents) {
const G4Material* mat = x.first;
G4int nn = mat->GetNumberOfElements();
for(G4int j=0; j<nn; ++j) {
G4bool yes = true;
const G4Element* elm = mat->GetElement(j);
for(G4int k=0; k<fNumberOfElements; ++k) {
if(elm == (*theElementVector)[k]) {
yes = false;
break;
}
}
if(yes) { ++nel; }
}
}
// resize vectors
if(nel > fNbComponents) {
delete [] fAtomsVector;
fAtomsVector = new G4int[nel];
G4double* v = new G4double[nel];
for(G4int i=0; i<fNumberOfElements; ++i) {
fAtomsVector[i] = 0;
v[i] = fMassFractionVector[i];
}
delete [] fMassFractionVector;
fMassFractionVector = v;
for(G4int i=fNumberOfElements; i<nel; ++i) {
fAtomsVector[i] = 0;
fMassFractionVector[i] = 0.0;
}
}
// filling
for(auto & x : fMatComponents) {
const G4Material* mat = x.first;
G4double frac = x.second;
G4int nn = mat->GetNumberOfElements();
const G4double* elmFrac = mat->GetFractionVector();
for(G4int j=0; j<nn; ++j) {
G4bool yes = true;
const G4Element* elm = mat->GetElement(j);
for(G4int k=0; k<fNumberOfElements; ++k) {
if(elm == (*theElementVector)[k]) {
fMassFractionVector[k] += frac*elmFrac[j];
yes = false;
break;
}
}
if(yes) {
theElementVector->push_back(const_cast<G4Element*>(elm));
fMassFractionVector[fNumberOfElements] = frac*elmFrac[j];
++fNumberOfElements;
}
}
}
}

// check sum of weights -- OK?
theElementVector = new G4ElementVector();
theElementVector->reserve(fNumberOfElements);
fAtomsVector = new G4int[fNumberOfElements];
fMassFractionVector = new G4double[fNumberOfElements];
G4double wtSum(0.0);
for (G4int i=0; i<fNumberOfElements; ++i) {
theElementVector->push_back((*fElm)[i]);
fMassFractionVector[i] = (*fElmFrac)[i];
wtSum += fMassFractionVector[i];
}
delete fElmFrac;
delete fElm;

// check sum of weights -- OK?
if (std::abs(1.-wtSum) > perThousand) {
G4ExceptionDescription ed;
ed << "For material " << fName << " sum of fractional masses "
Expand Down
2 changes: 1 addition & 1 deletion source/materials/src/G4NistMaterialBuilder.cc
View file
Expand Up @@ -734,7 +734,7 @@ void G4NistMaterialBuilder::NistSimpleMaterials()
AddMaterial("G4_Li", 0.534 , 3, 40. );
AddMaterial("G4_Be", 1.848 , 4, 63.7);
AddMaterial("G4_B" , 2.37 , 5, 76. );
AddMaterial("G4_C" , 2. , 6, 81. );
AddMaterial("G4_C" , 2. , 6, 78. );
AddMaterial("G4_N" , 1.16520e-3, 7, 82. , 1, kStateGas);
AddMaterial("G4_O" , 1.33151e-3, 8, 95. , 1, kStateGas);
AddMaterial("G4_F" , 1.58029e-3, 9, 115. , 1, kStateGas);
Expand Down

0 comments on commit 48e04df

Please sign in to comment.