Geant4 Cross Reference |
1 // 1 // 2 // ******************************************* 2 // ******************************************************************** 3 // * License and Disclaimer 3 // * License and Disclaimer * 4 // * 4 // * * 5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of * 6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and * 7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file * 8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These * 9 // * include a list of copyright holders. 9 // * include a list of copyright holders. * 10 // * 10 // * * 11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing * 12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this * 13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, * 14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its * 15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above * 16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. * 17 // * 17 // * * 18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and * 19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. * 20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or * 21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its * 22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your * 23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. * 24 // ******************************************* 24 // ******************************************************************** 25 << 25 // >> 26 // $Id: G4Material.cc 94016 2015-11-05 10:14:49Z gcosmo $ >> 27 // >> 28 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... >> 29 // 26 // 26-06-96, Code uses operators (+=, *=, ++, 30 // 26-06-96, Code uses operators (+=, *=, ++, -> etc.) correctly, P. Urban 27 // 10-07-96, new data members added by L.Urban 31 // 10-07-96, new data members added by L.Urban 28 // 12-12-96, new data memberfFreeElecDensitys << 32 // 12-12-96, new data members added by L.Urban 29 // 20-01-97, aesthetic rearrangement. RadLengt 33 // 20-01-97, aesthetic rearrangement. RadLength calculation modified. 30 // Data members Zeff and Aeff REMOVE 34 // Data members Zeff and Aeff REMOVED (i.e. passed to the Elements). 31 // (local definition of Zeff in Dens 35 // (local definition of Zeff in DensityEffect and FluctModel...) 32 // Vacuum defined as a G4State. Mixt << 36 // Vacuum defined as a G4State. Mixture flag removed, M.Maire. 33 // 29-01-97, State=Vacuum automatically set de 37 // 29-01-97, State=Vacuum automatically set density=0 in the contructors. 34 // Subsequent protections have been << 38 // Subsequent protections have been put in the calculation of 35 // MeanExcEnergy, ShellCorrectionVec 39 // MeanExcEnergy, ShellCorrectionVector, DensityEffect, M.Maire. 36 // 11-02-97, ComputeDensityEffect() rearranged 40 // 11-02-97, ComputeDensityEffect() rearranged, M.Maire. 37 // 20-03-97, corrected initialization of point 41 // 20-03-97, corrected initialization of pointers, M.Maire. 38 // 28-05-98, the kState=kVacuum has been remov 42 // 28-05-98, the kState=kVacuum has been removed. 39 // automatic check for a minimal den << 43 // automatic check for a minimal density, M.Maire 40 // 12-06-98, new method AddMaterial() allowing << 44 // 12-06-98, new method AddMaterial() allowing mixture of materials, M.Maire 41 // 09-07-98, ionisation parameters removed fro 45 // 09-07-98, ionisation parameters removed from the class, M.Maire 42 // 05-10-98, change names: NumDensity -> NbOfA 46 // 05-10-98, change names: NumDensity -> NbOfAtomsPerVolume 43 // 18-11-98, new interface to SandiaTable 47 // 18-11-98, new interface to SandiaTable 44 // 19-01-99 enlarge tolerance on test of cohe 48 // 19-01-99 enlarge tolerance on test of coherence of gas conditions 45 // 19-07-99, Constructors with chemicalFormula 49 // 19-07-99, Constructors with chemicalFormula added by V.Ivanchenko 46 // 16-01-01, Nuclear interaction length, M.Mai 50 // 16-01-01, Nuclear interaction length, M.Maire 47 // 12-03-01, G4bool fImplicitElement; 51 // 12-03-01, G4bool fImplicitElement; 48 // copy constructor and assignement 52 // copy constructor and assignement operator revised (mma) 49 // 03-05-01, flux.precision(prec) at begin/end 53 // 03-05-01, flux.precision(prec) at begin/end of operator<< 50 // 17-07-01, migration to STL. M. Verderi. 54 // 17-07-01, migration to STL. M. Verderi. 51 // 14-09-01, Suppression of the data member fI 55 // 14-09-01, Suppression of the data member fIndexInTable 52 // 26-02-02, fIndexInTable renewed 56 // 26-02-02, fIndexInTable renewed 53 // 16-04-02, G4Exception put in constructor wi 57 // 16-04-02, G4Exception put in constructor with chemical formula 54 // 06-05-02, remove the check of the ideal gas 58 // 06-05-02, remove the check of the ideal gas state equation 55 // 06-08-02, remove constructors with chemical 59 // 06-08-02, remove constructors with chemical formula (mma) 56 // 22-01-04, proper STL handling of theElement 60 // 22-01-04, proper STL handling of theElementVector (Hisaya) 57 // 30-03-05, warning in GetMaterial(materialNa << 61 // 30-03-05, warning in GetMaterial(materialName) 58 // 09-03-06, minor change of printout (V.Ivanc << 62 // 09-03-06, minor change of printout (V.Ivanchenko) 59 // 10-01-07, compute fAtomVector in the case o << 63 // 10-01-07, compute fAtomVector in the case of mass fraction (V.Ivanchenko) 60 // 27-07-07, improve destructor (V.Ivanchenko) << 64 // 27-07-07, improve destructor (V.Ivanchenko) 61 // 18-10-07, moved definition of mat index to << 65 // 18-10-07, move definition of material index to InitialisePointers (V.Ivanchenko) 62 // 13-08-08, do not use fixed size arrays (V.I 66 // 13-08-08, do not use fixed size arrays (V.Ivanchenko) 63 // 26-10-11, new scheme for G4Exception (mma) 67 // 26-10-11, new scheme for G4Exception (mma) 64 // 13-04-12, map<G4Material*,G4double> fMatCom 68 // 13-04-12, map<G4Material*,G4double> fMatComponents, filled in AddMaterial() 65 // 21-04-12, fMassOfMolecule, computed for Ato 69 // 21-04-12, fMassOfMolecule, computed for AtomsCount (mma) >> 70 // >> 71 // >> 72 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 66 73 67 #include "G4Material.hh" << 74 #include <iomanip> 68 75 69 #include "G4ApplicationState.hh" << 76 #include "G4Material.hh" 70 #include "G4AtomicShells.hh" << 71 #include "G4ExtendedMaterial.hh" << 72 #include "G4Pow.hh" << 73 #include "G4NistManager.hh" 77 #include "G4NistManager.hh" >> 78 #include "G4UnitsTable.hh" 74 #include "G4PhysicalConstants.hh" 79 #include "G4PhysicalConstants.hh" 75 #include "G4StateManager.hh" << 76 #include "G4SystemOfUnits.hh" 80 #include "G4SystemOfUnits.hh" 77 #include "G4UnitsTable.hh" << 81 #include "G4Exp.hh" >> 82 #include "G4Log.hh" 78 83 79 #include <iomanip> << 84 G4MaterialTable G4Material::theMaterialTable; 80 85 81 //....oooOO0OOooo........oooOO0OOooo........oo 86 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 82 87 83 // Constructor to create a material from scrat 88 // Constructor to create a material from scratch 84 89 85 G4Material::G4Material(const G4String& name, G << 90 G4Material::G4Material(const G4String& name, G4double z, 86 G4State state, G4double temp, G4double press << 91 G4double a, G4double density, 87 : fName(name) << 92 G4State state, G4double temp, G4double pressure) >> 93 : fName(name) 88 { 94 { 89 InitializePointers(); 95 InitializePointers(); 90 << 96 91 if (density < universe_mean_density) { << 97 if (density < universe_mean_density) 92 G4cout << " G4Material WARNING:" << 98 { 93 << " define a material with density << 99 G4cout << " G4Material WARNING:" 94 << " The material " << name << " wi << 100 << " define a material with density=0 is not allowed. \n" 95 << " default minimal density: " << << 101 << " The material " << name << " will be constructed with the" 96 << G4endl; << 102 << " default minimal density: " << universe_mean_density/(g/cm3) 97 density = universe_mean_density; << 103 << "g/cm3" << G4endl; 98 } << 104 density = universe_mean_density; 99 << 105 } 100 fDensity = density; << 106 101 fState = state; << 107 fDensity = density; 102 fTemp = temp; << 108 fState = state; >> 109 fTemp = temp; 103 fPressure = pressure; 110 fPressure = pressure; 104 111 105 // Initialize theElementVector allocating on 112 // Initialize theElementVector allocating one 106 // element corresponding to this material 113 // element corresponding to this material 107 fNbComponents = fNumberOfElements = 1; << 114 maxNbComponents = fNumberOfComponents = fNumberOfElements = 1; 108 theElementVector = new G4ElementVector(); << 115 fArrayLength = maxNbComponents; 109 << 116 theElementVector = new G4ElementVector(); 110 // take element from DB << 117 111 G4NistManager* nist = G4NistManager::Instanc << 118 const std::vector<G4String> elmnames = >> 119 G4NistManager::Instance()->GetNistElementNames(); >> 120 G4String enam, snam; 112 G4int iz = G4lrint(z); 121 G4int iz = G4lrint(z); 113 auto elm = nist->FindOrBuildElement(iz); << 122 if(iz < (G4int)elmnames.size()) { 114 if (elm == nullptr) { << 123 snam = elmnames[iz]; 115 elm = new G4Element("ELM_" + name, name, z << 124 enam = snam; 116 } << 125 } else { 117 theElementVector->push_back(elm); << 126 enam = "ELM_" + name; 118 << 127 snam = name; 119 fMassFractionVector = new G4double[1]; << 128 } 120 fMassFractionVector[0] = 1.; << 129 theElementVector->push_back(new G4Element(enam, snam, z, a)); 121 fMassOfMolecule = a / CLHEP::Avogadro; << 130 122 << 131 fMassFractionVector = new G4double[1]; 123 if (fState == kStateUndefined) { << 132 fMassFractionVector[0] = 1. ; 124 if (fDensity > kGasThreshold) { << 133 fMassOfMolecule = a/Avogadro; 125 fState = kStateSolid; << 134 126 } << 135 if (fState == kStateUndefined) 127 else { << 136 { 128 fState = kStateGas; << 137 if (fDensity > kGasThreshold) { fState = kStateSolid; } >> 138 else { fState = kStateGas; } 129 } 139 } 130 } << 131 140 132 ComputeDerivedQuantities(); 141 ComputeDerivedQuantities(); 133 } 142 } 134 143 135 //....oooOO0OOooo........oooOO0OOooo........oo 144 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 136 145 137 // Constructor to create a material from a Lis 146 // Constructor to create a material from a List of constituents 138 // (elements and/or materials) added with Add 147 // (elements and/or materials) added with AddElement or AddMaterial 139 148 140 G4Material::G4Material(const G4String& name, G << 149 G4Material::G4Material(const G4String& name, G4double density, 141 G4double temp, G4double pressure) << 150 G4int nComponents, 142 : fName(name) << 151 G4State state, G4double temp, G4double pressure) >> 152 : fName(name) 143 { 153 { 144 InitializePointers(); 154 InitializePointers(); 145 << 155 146 if (density < universe_mean_density) { << 156 if (density < universe_mean_density) 147 G4cout << "--- Warning from G4Material::G4 << 157 { 148 << " define a material with density << 158 G4cout << "--- Warning from G4Material::G4Material()" 149 << " The material " << name << " wi << 159 << " define a material with density=0 is not allowed. \n" 150 << " default minimal density: " << << 160 << " The material " << name << " will be constructed with the" 151 << G4endl; << 161 << " default minimal density: " << universe_mean_density/(g/cm3) 152 density = universe_mean_density; << 162 << "g/cm3" << G4endl; 153 } << 163 density = universe_mean_density; 154 << 164 } 155 fDensity = density; << 165 156 fState = state; << 166 fDensity = density; 157 fTemp = temp; << 167 fState = state; >> 168 fTemp = temp; 158 fPressure = pressure; 169 fPressure = pressure; 159 << 170 160 fNbComponents = nComponents; << 171 maxNbComponents = nComponents; 161 fMassFraction = true; << 172 fArrayLength = maxNbComponents; 162 << 173 fNumberOfComponents = fNumberOfElements = 0; 163 if (fState == kStateUndefined) { << 174 theElementVector = new G4ElementVector(); 164 if (fDensity > kGasThreshold) { << 175 theElementVector->reserve(maxNbComponents); 165 fState = kStateSolid; << 176 166 } << 177 if (fState == kStateUndefined) 167 else { << 178 { 168 fState = kStateGas; << 179 if (fDensity > kGasThreshold) { fState = kStateSolid; } >> 180 else { fState = kStateGas; } 169 } 181 } 170 } << 171 } 182 } 172 183 173 //....oooOO0OOooo........oooOO0OOooo........oo 184 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 174 185 175 // Constructor to create a material from base 186 // Constructor to create a material from base material 176 187 177 G4Material::G4Material(const G4String& name, G << 188 G4Material::G4Material(const G4String& name, G4double density, 178 G4State state, G4double temp, G4double press << 189 const G4Material* bmat, 179 : fName(name) << 190 G4State state, G4double temp, G4double pressure) >> 191 : fName(name) 180 { 192 { 181 InitializePointers(); 193 InitializePointers(); >> 194 >> 195 if (density < universe_mean_density) >> 196 { >> 197 G4cout << "--- Warning from G4Material::G4Material()" >> 198 << " define a material with density=0 is not allowed. \n" >> 199 << " The material " << name << " will be constructed with the" >> 200 << " default minimal density: " << universe_mean_density/(g/cm3) >> 201 << "g/cm3" << G4endl; >> 202 density = universe_mean_density; >> 203 } 182 204 183 if (density < universe_mean_density) { << 205 fDensity = density; 184 G4cout << "--- Warning from G4Material::G4 << 206 fState = state; 185 << " define a material with density << 207 fTemp = temp; 186 << " The material " << name << " wi << 187 << " default minimal density: " << << 188 << G4endl; << 189 density = universe_mean_density; << 190 } << 191 << 192 fDensity = density; << 193 fState = state; << 194 fTemp = temp; << 195 fPressure = pressure; 208 fPressure = pressure; 196 209 197 fBaseMaterial = bmat; 210 fBaseMaterial = bmat; 198 auto ptr = bmat; << 199 if (nullptr != ptr) { << 200 while (true) { << 201 ptr = ptr->GetBaseMaterial(); << 202 if (nullptr == ptr) { << 203 break; << 204 } << 205 fBaseMaterial = ptr; << 206 } << 207 } << 208 << 209 fChemicalFormula = fBaseMaterial->GetChemica 211 fChemicalFormula = fBaseMaterial->GetChemicalFormula(); 210 fMassOfMolecule = fBaseMaterial->GetMassOfMo << 212 fMassOfMolecule = fBaseMaterial->GetMassOfMolecule(); 211 213 212 fNumberOfElements = (G4int)fBaseMaterial->Ge << 214 fNumberOfElements = fBaseMaterial->GetNumberOfElements(); 213 fNbComponents = fNumberOfElements; << 215 maxNbComponents = fNumberOfElements; >> 216 fNumberOfComponents = fNumberOfElements; >> 217 >> 218 fMaterialPropertiesTable = fBaseMaterial->GetMaterialPropertiesTable(); 214 219 215 CopyPointersOfBaseMaterial(); 220 CopyPointersOfBaseMaterial(); 216 } 221 } 217 222 218 //....oooOO0OOooo........oooOO0OOooo........oo 223 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 219 224 >> 225 // Fake default constructor - sets only member data and allocates memory >> 226 // for usage restricted to object persistency >> 227 >> 228 G4Material::G4Material(__void__&) >> 229 : fChemicalFormula(""), fDensity(0.0), fState(kStateUndefined), fTemp(0.0), >> 230 fPressure(0.0), maxNbComponents(0), fArrayLength(0), >> 231 fNumberOfComponents(0), fNumberOfElements(0), theElementVector(0), >> 232 fMassFractionVector(0), fAtomsVector(0), fMaterialPropertiesTable(0), >> 233 fIndexInTable(0), VecNbOfAtomsPerVolume(0), TotNbOfAtomsPerVolume(0), >> 234 TotNbOfElectPerVolume(0), fRadlen(0.0), fNuclInterLen(0.0), fIonisation(0), >> 235 fSandiaTable(0), fBaseMaterial(0), fMassOfMolecule(0.0) >> 236 { >> 237 } >> 238 >> 239 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... >> 240 220 G4Material::~G4Material() 241 G4Material::~G4Material() 221 { 242 { 222 if (fBaseMaterial == nullptr) { << 243 // G4cout << "### Destruction of material " << fName << " started" <<G4endl; 223 delete theElementVector; << 244 if(!fBaseMaterial) { 224 delete fSandiaTable; << 245 if (theElementVector) { delete theElementVector; } 225 delete[] fMassFractionVector; << 246 if (fMassFractionVector) { delete [] fMassFractionVector; } 226 delete[] fAtomsVector; << 247 if (fAtomsVector) { delete [] fAtomsVector; } >> 248 if (fSandiaTable) { delete fSandiaTable; } 227 } 249 } 228 delete fIonisation; << 250 if (fIonisation) { delete fIonisation; } 229 delete[] fVecNbOfAtomsPerVolume; << 251 if (VecNbOfAtomsPerVolume) { delete [] VecNbOfAtomsPerVolume; } 230 252 231 // Remove this material from the MaterialTab << 253 // Remove this material from theMaterialTable. 232 // 254 // 233 (*GetMaterialTable())[fIndexInTable] = nullp << 255 theMaterialTable[fIndexInTable] = 0; 234 } 256 } 235 257 236 //....oooOO0OOooo........oooOO0OOooo........oo 258 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 237 259 238 void G4Material::InitializePointers() 260 void G4Material::InitializePointers() 239 { 261 { 240 fBaseMaterial = nullptr; << 262 theElementVector = 0; 241 fMaterialPropertiesTable = nullptr; << 263 fMassFractionVector = 0; 242 theElementVector = nullptr; << 264 fAtomsVector = 0; 243 fAtomsVector = nullptr; << 265 fMaterialPropertiesTable = 0; 244 fMassFractionVector = nullptr; << 266 245 fVecNbOfAtomsPerVolume = nullptr; << 267 VecNbOfAtomsPerVolume = 0; 246 << 268 fBaseMaterial = 0; 247 fIonisation = nullptr; << 269 248 fSandiaTable = nullptr; << 270 fChemicalFormula = ""; 249 << 271 250 fDensity = fFreeElecDensity = fTemp = fPress << 272 // initilized data members 251 fTotNbOfAtomsPerVolume = 0.0; << 273 fDensity = 0.0; 252 fTotNbOfElectPerVolume = 0.0; << 274 fState = kStateUndefined; 253 fRadlen = fNuclInterLen = fMassOfMolecule = << 275 fTemp = 0.0; 254 << 276 fPressure = 0.0; 255 fState = kStateUndefined; << 277 maxNbComponents = 0; 256 fNumberOfElements = fNbComponents = fIdxComp << 278 fArrayLength = 0; 257 fMassFraction = true; << 279 TotNbOfAtomsPerVolume = 0; 258 fChemicalFormula = ""; << 280 TotNbOfElectPerVolume = 0; >> 281 fRadlen = 0.0; >> 282 fNuclInterLen = 0.0; >> 283 fMassOfMolecule = 0.0; >> 284 >> 285 fIonisation = 0; >> 286 fSandiaTable = 0; 259 287 260 // Store in the static Table of Materials 288 // Store in the static Table of Materials 261 fIndexInTable = GetMaterialTable()->size(); << 289 fIndexInTable = theMaterialTable.size(); 262 for (std::size_t i = 0; i < fIndexInTable; + << 290 for(size_t i=0; i<fIndexInTable; ++i) { 263 if ((*GetMaterialTable())[i]->GetName() == << 291 if(theMaterialTable[i]->GetName() == fName) { 264 G4cout << "G4Material WARNING: duplicate << 292 G4cout << "G4Material WARNING: duplicate name of material " >> 293 << fName << G4endl; 265 break; 294 break; 266 } 295 } 267 } 296 } 268 GetMaterialTable()->push_back(this); << 297 theMaterialTable.push_back(this); 269 } 298 } 270 299 271 //....oooOO0OOooo........oooOO0OOooo........oo 300 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 272 301 273 void G4Material::ComputeDerivedQuantities() 302 void G4Material::ComputeDerivedQuantities() 274 { 303 { 275 // Header routine to compute various propert 304 // Header routine to compute various properties of material. 276 // << 305 // >> 306 277 // Number of atoms per volume (per element), 307 // Number of atoms per volume (per element), total nb of electrons per volume 278 G4double Zi, Ai; 308 G4double Zi, Ai; 279 fTotNbOfAtomsPerVolume = 0.; << 309 TotNbOfAtomsPerVolume = 0.; 280 delete[] fVecNbOfAtomsPerVolume; << 310 if (VecNbOfAtomsPerVolume) { delete [] VecNbOfAtomsPerVolume; } 281 fVecNbOfAtomsPerVolume = new G4double[fNumbe << 311 VecNbOfAtomsPerVolume = new G4double[fNumberOfElements]; 282 fTotNbOfElectPerVolume = 0.; << 312 TotNbOfElectPerVolume = 0.; 283 fFreeElecDensity = 0.; << 313 for (G4int i=0; i<fNumberOfElements; ++i) { 284 const G4double elecTh = 15. * CLHEP::eV; // << 314 Zi = (*theElementVector)[i]->GetZ(); 285 for (G4int i = 0; i < fNumberOfElements; ++i << 315 Ai = (*theElementVector)[i]->GetA(); 286 Zi = (*theElementVector)[i]->GetZ(); << 316 VecNbOfAtomsPerVolume[i] = Avogadro*fDensity*fMassFractionVector[i]/Ai; 287 Ai = (*theElementVector)[i]->GetA(); << 317 TotNbOfAtomsPerVolume += VecNbOfAtomsPerVolume[i]; 288 fVecNbOfAtomsPerVolume[i] = Avogadro * fDe << 318 TotNbOfElectPerVolume += VecNbOfAtomsPerVolume[i]*Zi; 289 fTotNbOfAtomsPerVolume += fVecNbOfAtomsPer << 290 fTotNbOfElectPerVolume += fVecNbOfAtomsPer << 291 if (fState != kStateGas) { << 292 fFreeElecDensity += << 293 fVecNbOfAtomsPerVolume[i] * G4AtomicSh << 294 } << 295 } 319 } 296 << 320 297 ComputeRadiationLength(); 321 ComputeRadiationLength(); 298 ComputeNuclearInterLength(); 322 ComputeNuclearInterLength(); 299 323 300 if (fIonisation == nullptr) { << 324 if (fIonisation) { delete fIonisation; } 301 fIonisation = new G4IonisParamMat(this); << 325 fIonisation = new G4IonisParamMat(this); 302 } << 326 if (fSandiaTable) { delete fSandiaTable; } 303 if (fSandiaTable == nullptr) { << 327 fSandiaTable = new G4SandiaTable(this); 304 fSandiaTable = new G4SandiaTable(this); << 305 } << 306 } 328 } 307 329 308 //....oooOO0OOooo........oooOO0OOooo........oo 330 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 309 331 310 void G4Material::CopyPointersOfBaseMaterial() 332 void G4Material::CopyPointersOfBaseMaterial() 311 { 333 { 312 G4double factor = fDensity / fBaseMaterial-> << 334 G4double factor = fDensity/fBaseMaterial->GetDensity(); 313 fTotNbOfAtomsPerVolume = factor * fBaseMater << 335 TotNbOfAtomsPerVolume = factor*fBaseMaterial->GetTotNbOfAtomsPerVolume(); 314 fTotNbOfElectPerVolume = factor * fBaseMater << 336 TotNbOfElectPerVolume = factor*fBaseMaterial->GetTotNbOfElectPerVolume(); 315 fFreeElecDensity = factor * fBaseMaterial->G << 337 316 << 338 theElementVector = 317 if (fState == kStateUndefined) { << 339 const_cast<G4ElementVector*>(fBaseMaterial->GetElementVector()); 318 fState = fBaseMaterial->GetState(); << 340 fMassFractionVector = 319 } << 341 const_cast<G4double*>(fBaseMaterial->GetFractionVector()); 320 << 321 theElementVector = const_cast<G4ElementVecto << 322 fMassFractionVector = const_cast<G4double*>( << 323 fAtomsVector = const_cast<G4int*>(fBaseMater 342 fAtomsVector = const_cast<G4int*>(fBaseMaterial->GetAtomsVector()); 324 343 325 const G4double* v = fBaseMaterial->GetVecNbO 344 const G4double* v = fBaseMaterial->GetVecNbOfAtomsPerVolume(); 326 delete[] fVecNbOfAtomsPerVolume; << 345 if (VecNbOfAtomsPerVolume) { delete [] VecNbOfAtomsPerVolume; } 327 fVecNbOfAtomsPerVolume = new G4double[fNumbe << 346 VecNbOfAtomsPerVolume = new G4double[fNumberOfElements]; 328 for (G4int i = 0; i < fNumberOfElements; ++i << 347 for (G4int i=0; i<fNumberOfElements; ++i) { 329 fVecNbOfAtomsPerVolume[i] = factor * v[i]; << 348 VecNbOfAtomsPerVolume[i] = factor*v[i]; 330 } 349 } 331 fRadlen = fBaseMaterial->GetRadlen() / facto << 350 fRadlen = fBaseMaterial->GetRadlen()/factor; 332 fNuclInterLen = fBaseMaterial->GetNuclearInt << 351 fNuclInterLen = fBaseMaterial->GetNuclearInterLength()/factor; 333 352 334 if (fIonisation == nullptr) { << 353 if (fIonisation) { delete fIonisation; } 335 fIonisation = new G4IonisParamMat(this); << 354 fIonisation = new G4IonisParamMat(this); 336 } << 337 fIonisation->SetMeanExcitationEnergy(fBaseMa << 338 if (fBaseMaterial->GetIonisation()->GetDensi << 339 ComputeDensityEffectOnFly(true); << 340 } << 341 355 342 fSandiaTable = fBaseMaterial->GetSandiaTable 356 fSandiaTable = fBaseMaterial->GetSandiaTable(); 343 fMaterialPropertiesTable = fBaseMaterial->Ge << 357 fIonisation->SetMeanExcitationEnergy(fBaseMaterial->GetIonisation()->GetMeanExcitationEnergy()); 344 } 358 } 345 359 346 //....oooOO0OOooo........oooOO0OOooo........oo 360 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 347 361 348 void G4Material::AddElementByNumberOfAtoms(con << 362 // AddElement -- composition by atom count 349 { << 350 // perform checks consistency << 351 if (0 == fIdxComponent) { << 352 fMassFraction = false; << 353 fAtoms = new std::vector<G4int>; << 354 fElm = new std::vector<const G4Element*>; << 355 } << 356 if (fIdxComponent >= fNbComponents) { << 357 G4ExceptionDescription ed; << 358 ed << "For material " << fName << " and ad << 359 << " with Natoms=" << nAtoms << 360 << " wrong attempt to add more than the << 361 << " >= " << fNbComponents; << 362 G4Exception("G4Material::AddElementByNumbe << 363 } << 364 if (fMassFraction) { << 365 G4ExceptionDescription ed; << 366 ed << "For material " << fName << " and ad << 367 << " with Natoms=" << nAtoms << " probl << 368 << "addition of elements by mass fracti << 369 G4Exception("G4Material::AddElementByNumbe << 370 } << 371 if (0 >= nAtoms) { << 372 G4ExceptionDescription ed; << 373 ed << "For material " << fName << " and ad << 374 << " with Natoms=" << nAtoms << " probl << 375 G4Exception("G4Material::AddElementByNumbe << 376 } << 377 << 378 // filling << 379 G4bool isAdded = false; << 380 if (! fElm->empty()) { << 381 for (G4int i = 0; i < fNumberOfElements; + << 382 if (elm == (*fElm)[i]) { << 383 (*fAtoms)[i] += nAtoms; << 384 isAdded = true; << 385 break; << 386 } << 387 } << 388 } << 389 if (! isAdded) { << 390 fElm->push_back(elm); << 391 fAtoms->push_back(nAtoms); << 392 ++fNumberOfElements; << 393 } << 394 ++fIdxComponent; << 395 << 396 // is filled - complete composition of atoms << 397 if (fIdxComponent == fNbComponents) { << 398 theElementVector = new G4ElementVector(); << 399 theElementVector->reserve(fNumberOfElement << 400 fAtomsVector = new G4int[fNumberOfElements << 401 fMassFractionVector = new G4double[fNumber << 402 363 >> 364 void G4Material::AddElement(G4Element* element, G4int nAtoms) >> 365 { >> 366 // initialization >> 367 if ( fNumberOfElements == 0 ) { >> 368 fAtomsVector = new G4int [fArrayLength]; >> 369 fMassFractionVector = new G4double[fArrayLength]; >> 370 } >> 371 >> 372 // filling ... >> 373 if ( fNumberOfElements < maxNbComponents ) { >> 374 theElementVector->push_back(element); >> 375 fAtomsVector[fNumberOfElements] = nAtoms; >> 376 fNumberOfComponents = ++fNumberOfElements; >> 377 } else { >> 378 G4cout << "G4Material::AddElement ERROR for " << fName << " nElement= " >> 379 << fNumberOfElements << G4endl; >> 380 G4Exception ("G4Material::AddElement()", "mat031", FatalException, >> 381 "Attempt to add more than the declared number of elements."); >> 382 } >> 383 // filled. >> 384 if ( fNumberOfElements == maxNbComponents ) { >> 385 // compute proportion by mass >> 386 G4int i=0; 403 G4double Amol = 0.; 387 G4double Amol = 0.; 404 for (G4int i = 0; i < fNumberOfElements; + << 388 for (i=0; i<fNumberOfElements; ++i) { 405 theElementVector->push_back((*fElm)[i]); << 389 G4double w = fAtomsVector[i]*(*theElementVector)[i]->GetA(); 406 fAtomsVector[i] = (*fAtoms)[i]; << 407 G4double w = fAtomsVector[i] * (*fElm)[i << 408 Amol += w; 390 Amol += w; 409 fMassFractionVector[i] = w; 391 fMassFractionVector[i] = w; 410 } 392 } 411 for (G4int i = 0; i < fNumberOfElements; + << 393 for (i=0; i<fNumberOfElements; ++i) { 412 fMassFractionVector[i] /= Amol; 394 fMassFractionVector[i] /= Amol; 413 } 395 } 414 delete fAtoms; << 396 415 delete fElm; << 397 fMassOfMolecule = Amol/Avogadro; 416 fMassOfMolecule = Amol / CLHEP::Avogadro; << 417 ComputeDerivedQuantities(); 398 ComputeDerivedQuantities(); 418 } 399 } 419 } 400 } 420 401 421 //....oooOO0OOooo........oooOO0OOooo........oo 402 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 422 403 423 void G4Material::AddElementByMassFraction(cons << 404 // AddElement -- composition by fraction of mass 424 { << 425 // perform checks consistency << 426 if (fraction < 0.0 || fraction > 1.0) { << 427 G4ExceptionDescription ed; << 428 ed << "For material " << fName << " and ad << 429 << " massFraction= " << fraction << " i << 430 G4Exception("G4Material::AddElementByMassF << 431 } << 432 if (! fMassFraction) { << 433 G4ExceptionDescription ed; << 434 ed << "For material " << fName << " and ad << 435 << ", massFraction= " << fraction << ", << 436 << " problem: cannot add by mass fracti << 437 << "addition of elements by number of a << 438 G4Exception("G4Material::AddElementByMassF << 439 } << 440 if (fIdxComponent >= fNbComponents) { << 441 G4ExceptionDescription ed; << 442 ed << "For material " << fName << " and ad << 443 << ", massFraction= " << fraction << ", << 444 << "; attempt to add more than the decl << 445 << " >= " << fNbComponents; << 446 G4Exception("G4Material::AddElementByMassF << 447 } << 448 if (0 == fIdxComponent) { << 449 fElmFrac = new std::vector<G4double>; << 450 fElm = new std::vector<const G4Element*>; << 451 } << 452 405 453 // filling << 406 void G4Material::AddElement(G4Element* element, G4double fraction) 454 G4bool isAdded = false; << 407 { 455 if (! fElm->empty()) { << 408 if(fraction < 0.0 || fraction > 1.0) { 456 for (G4int i = 0; i < fNumberOfElements; + << 409 G4cout << "G4Material::AddElement ERROR for " << fName << " and " 457 if (elm == (*fElm)[i]) { << 410 << element->GetName() << " mass fraction= " << fraction 458 (*fElmFrac)[i] += fraction; << 411 << " is wrong " << G4endl; 459 isAdded = true; << 412 G4Exception ("G4Material::AddElement()", "mat032", FatalException, 460 break; << 413 "Attempt to add element with wrong mass fraction"); 461 } << 414 } >> 415 // initialization >> 416 if (fNumberOfComponents == 0) { >> 417 fMassFractionVector = new G4double[fArrayLength]; >> 418 fAtomsVector = new G4int [fArrayLength]; >> 419 } >> 420 // filling ... >> 421 if (fNumberOfComponents < maxNbComponents) { >> 422 G4int el = 0; >> 423 // Loop checking, 07-Aug-2015, Vladimir Ivanchenko >> 424 while ((el<fNumberOfElements)&&(element!=(*theElementVector)[el])) { ++el; } >> 425 if (el<fNumberOfElements) fMassFractionVector[el] += fraction; >> 426 else { >> 427 theElementVector->push_back(element); >> 428 fMassFractionVector[el] = fraction; >> 429 ++fNumberOfElements; 462 } 430 } 463 } << 431 ++fNumberOfComponents; 464 if (! isAdded) { << 432 } else { 465 fElm->push_back(elm); << 433 G4cout << "G4Material::AddElement ERROR for " << fName << " nElement= " 466 fElmFrac->push_back(fraction); << 434 << fNumberOfElements << G4endl; 467 ++fNumberOfElements; << 435 G4Exception ("G4Material::AddElement()", "mat033", FatalException, 468 } << 436 "Attempt to add more than the declared number of elements."); 469 ++fIdxComponent; << 437 } 470 << 438 471 // is filled << 439 // filled. 472 if (fIdxComponent == fNbComponents) { << 440 if (fNumberOfComponents == maxNbComponents) { 473 FillVectors(); << 441 >> 442 G4int i=0; >> 443 G4double Zmol(0.), Amol(0.); >> 444 // check sum of weights -- OK? >> 445 G4double wtSum(0.0); >> 446 for (i=0; i<fNumberOfElements; ++i) { >> 447 wtSum += fMassFractionVector[i]; >> 448 Zmol += fMassFractionVector[i]*(*theElementVector)[i]->GetZ(); >> 449 Amol += fMassFractionVector[i]*(*theElementVector)[i]->GetA(); >> 450 } >> 451 if (std::fabs(1.-wtSum) > perThousand) { >> 452 G4cerr << "WARNING !! for " << fName << " sum of fractional masses " >> 453 << wtSum << " is not 1 - results may be wrong" >> 454 << G4endl; >> 455 } >> 456 for (i=0; i<fNumberOfElements; ++i) { >> 457 fAtomsVector[i] = >> 458 G4lrint(fMassFractionVector[i]*Amol/(*theElementVector)[i]->GetA()); >> 459 } >> 460 >> 461 ComputeDerivedQuantities(); 474 } 462 } 475 } 463 } 476 464 477 //....oooOO0OOooo........oooOO0OOooo........oo 465 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 478 466 479 // composition by fraction of mass << 467 // AddMaterial -- composition by fraction of mass >> 468 480 void G4Material::AddMaterial(G4Material* mater 469 void G4Material::AddMaterial(G4Material* material, G4double fraction) 481 { 470 { 482 if (fraction < 0.0 || fraction > 1.0) { << 471 if(fraction < 0.0 || fraction > 1.0) { 483 G4ExceptionDescription ed; << 472 G4cout << "G4Material::AddMaterial ERROR for " << fName << " and " 484 ed << "For material " << fName << " and ad << 473 << material->GetName() << " mass fraction= " << fraction 485 << ", massFraction= " << fraction << " << 474 << " is wrong "; 486 G4Exception("G4Material::AddMaterial()", " << 475 G4Exception ("G4Material::AddMaterial()", "mat034", FatalException, 487 } << 476 "Attempt to add material with wrong mass fraction"); 488 if (! fMassFraction) { << 477 } 489 G4ExceptionDescription ed; << 478 // initialization 490 ed << "For material " << fName << " and ad << 479 if (fNumberOfComponents == 0) { 491 << ", massFraction= " << fraction << ", << 480 fMassFractionVector = new G4double[fArrayLength]; 492 << " problem: cannot add by mass fracti << 481 fAtomsVector = new G4int [fArrayLength]; 493 << "addition of elements by number of a << 482 } 494 G4Exception("G4Material::AddMaterial()", " << 483 495 } << 484 G4int nelm = material->GetNumberOfElements(); 496 if (fIdxComponent >= fNbComponents) { << 485 497 G4ExceptionDescription ed; << 486 // arrays should be extended 498 ed << "For material " << fName << " and ad << 487 if(nelm > 1) { 499 << ", massFraction= " << fraction << 488 G4int nold = fArrayLength; 500 << "; attempt to add more than the decl << 489 fArrayLength += nelm - 1; 501 << " >= " << fNbComponents; << 490 G4double* v1 = new G4double[fArrayLength]; 502 G4Exception("G4Material::AddMaterial()", " << 491 G4int* i1 = new G4int[fArrayLength]; 503 } << 492 for(G4int i=0; i<nold; ++i) { 504 if (0 == fIdxComponent) { << 493 v1[i] = fMassFractionVector[i]; 505 fElmFrac = new std::vector<G4double>; << 494 i1[i] = fAtomsVector[i]; 506 fElm = new std::vector<const G4Element*>; << 495 } 507 } << 496 delete [] fAtomsVector; 508 << 497 delete [] fMassFractionVector; 509 // filling << 498 fMassFractionVector = v1; 510 auto nelm = (G4int)material->GetNumberOfElem << 499 fAtomsVector = i1; 511 for (G4int j = 0; j < nelm; ++j) { << 500 } 512 auto elm = material->GetElement(j); << 501 513 auto frac = material->GetFractionVector(); << 502 // filling ... 514 G4bool isAdded = false; << 503 if (fNumberOfComponents < maxNbComponents) { 515 if (! fElm->empty()) { << 504 for (G4int elm=0; elm<nelm; ++elm) 516 for (G4int i = 0; i < fNumberOfElements; << 505 { 517 if (elm == (*fElm)[i]) { << 506 G4Element* element = (*(material->GetElementVector()))[elm]; 518 (*fElmFrac)[i] += fraction * frac[j] << 507 G4int el = 0; 519 isAdded = true; << 508 // Loop checking, 07-Aug-2015, Vladimir Ivanchenko 520 break; << 509 while ((el<fNumberOfElements)&&(element!=(*theElementVector)[el])) el++; >> 510 if (el < fNumberOfElements) fMassFractionVector[el] += fraction >> 511 *(material->GetFractionVector())[elm]; >> 512 else { >> 513 theElementVector->push_back(element); >> 514 fMassFractionVector[el] = fraction >> 515 *(material->GetFractionVector())[elm]; >> 516 ++fNumberOfElements; 521 } 517 } 522 } << 518 } >> 519 ++fNumberOfComponents; >> 520 ///store massFraction of material component >> 521 fMatComponents[material] = fraction; >> 522 >> 523 } else { >> 524 G4cout << "G4Material::AddMaterial ERROR for " << fName << " nElement= " >> 525 << fNumberOfElements << G4endl; >> 526 G4Exception ("G4Material::AddMaterial()", "mat035", FatalException, >> 527 "Attempt to add more than the declared number of components."); >> 528 } >> 529 >> 530 // filled. >> 531 if (fNumberOfComponents == maxNbComponents) { >> 532 G4int i=0; >> 533 G4double Zmol(0.), Amol(0.); >> 534 // check sum of weights -- OK? >> 535 G4double wtSum(0.0); >> 536 for (i=0; i<fNumberOfElements; ++i) { >> 537 wtSum += fMassFractionVector[i]; >> 538 Zmol += fMassFractionVector[i]*(*theElementVector)[i]->GetZ(); >> 539 Amol += fMassFractionVector[i]*(*theElementVector)[i]->GetA(); >> 540 } >> 541 if (std::fabs(1.-wtSum) > perThousand) { >> 542 G4cout << "G4Material::AddMaterial WARNING !! for " << fName >> 543 << " sum of fractional masses " >> 544 << wtSum << " is not 1 - results may be wrong" >> 545 << G4endl; >> 546 } >> 547 for (i=0; i<fNumberOfElements; ++i) { >> 548 fAtomsVector[i] = >> 549 G4lrint(fMassFractionVector[i]*Amol/(*theElementVector)[i]->GetA()); 523 } 550 } 524 if (! isAdded) { << 551 525 fElm->push_back(elm); << 552 ComputeDerivedQuantities(); 526 fElmFrac->push_back(fraction * frac[j]); << 527 ++fNumberOfElements; << 528 } << 529 } << 530 << 531 fMatComponents[material] = fraction; << 532 ++fIdxComponent; << 533 << 534 // is filled << 535 if (fIdxComponent == fNbComponents) { << 536 FillVectors(); << 537 } << 538 } << 539 << 540 //....oooOO0OOooo........oooOO0OOooo........oo << 541 << 542 void G4Material::FillVectors() << 543 { << 544 // there are material components << 545 theElementVector = new G4ElementVector(); << 546 theElementVector->reserve(fNumberOfElements) << 547 fAtomsVector = new G4int[fNumberOfElements]; << 548 fMassFractionVector = new G4double[fNumberOf << 549 << 550 G4double wtSum(0.0); << 551 for (G4int i = 0; i < fNumberOfElements; ++i << 552 theElementVector->push_back((*fElm)[i]); << 553 fMassFractionVector[i] = (*fElmFrac)[i]; << 554 wtSum += fMassFractionVector[i]; << 555 } << 556 delete fElmFrac; << 557 delete fElm; << 558 << 559 // check sum of weights -- OK? << 560 if (std::abs(1. - wtSum) > perThousand) { << 561 G4ExceptionDescription ed; << 562 ed << "For material " << fName << " sum of << 563 << " is not 1 - results may be wrong"; << 564 G4Exception("G4Material::FillVectors()", " << 565 } << 566 G4double coeff = (wtSum > 0.0) ? 1. / wtSum << 567 G4double Amol(0.); << 568 for (G4int i = 0; i < fNumberOfElements; ++i << 569 fMassFractionVector[i] *= coeff; << 570 Amol += fMassFractionVector[i] * (*theElem << 571 } << 572 for (G4int i = 0; i < fNumberOfElements; ++i << 573 fAtomsVector[i] = G4lrint(fMassFractionVec << 574 } 553 } 575 ComputeDerivedQuantities(); << 576 } 554 } 577 555 578 //....oooOO0OOooo........oooOO0OOooo........oo 556 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 579 557 580 void G4Material::ComputeRadiationLength() 558 void G4Material::ComputeRadiationLength() 581 { 559 { 582 G4double radinv = 0.0; << 560 G4double radinv = 0.0 ; 583 for (G4int i = 0; i < fNumberOfElements; ++i << 561 for (G4int i=0;i<fNumberOfElements;++i) { 584 radinv += fVecNbOfAtomsPerVolume[i] * ((*t << 562 radinv += VecNbOfAtomsPerVolume[i]*((*theElementVector)[i]->GetfRadTsai()); 585 } 563 } 586 fRadlen = (radinv <= 0.0 ? DBL_MAX : 1. / ra << 564 fRadlen = (radinv <= 0.0 ? DBL_MAX : 1./radinv); 587 } 565 } 588 566 589 //....oooOO0OOooo........oooOO0OOooo........oo 567 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 590 568 591 void G4Material::ComputeNuclearInterLength() 569 void G4Material::ComputeNuclearInterLength() 592 { 570 { 593 const G4double lambda1 = CLHEP::amu*CLHEP::c << 571 static const G4double lambda0 = 35*CLHEP::g/CLHEP::cm2; >> 572 static const G4double twothird = 2.0/3.0; 594 G4double NILinv = 0.0; 573 G4double NILinv = 0.0; 595 for (G4int i = 0; i < fNumberOfElements; ++i << 574 for (G4int i=0; i<fNumberOfElements; ++i) { 596 G4int Z = (*theElementVector)[i]->GetZasIn << 575 G4int Z = G4lrint( (*theElementVector)[i]->GetZ()); 597 G4double A = (*theElementVector)[i]->GetN( 576 G4double A = (*theElementVector)[i]->GetN(); 598 if (1 == Z) { << 577 if(1 == Z) { 599 NILinv += fVecNbOfAtomsPerVolume[i] * A; << 578 NILinv += VecNbOfAtomsPerVolume[i]*A; 600 } << 579 } else { 601 else { << 580 NILinv += VecNbOfAtomsPerVolume[i]*G4Exp(twothird*G4Log(A)); 602 NILinv += fVecNbOfAtomsPerVolume[i] * G4 << 581 } 603 } << 604 } << 605 NILinv *= lambda1; << 606 fNuclInterLen = 1.e+20*CLHEP::m; << 607 if (fNuclInterLen * NILinv > 1.0) { << 608 fNuclInterLen = 1.0 / NILinv; << 609 } 582 } >> 583 NILinv *= amu/lambda0; >> 584 fNuclInterLen = (NILinv <= 0.0 ? DBL_MAX : 1./NILinv); 610 } 585 } 611 586 612 //....oooOO0OOooo........oooOO0OOooo........oo 587 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 613 588 614 void G4Material::SetChemicalFormula(const G4St << 589 G4MaterialTable* G4Material::GetMaterialTable() 615 { << 616 if (! IsLocked()) { << 617 fChemicalFormula = chF; << 618 } << 619 } << 620 << 621 //....oooOO0OOooo........oooOO0OOooo........oo << 622 << 623 void G4Material::SetFreeElectronDensity(G4doub << 624 { 590 { 625 if (val >= 0. && ! IsLocked()) { << 591 return &theMaterialTable; 626 fFreeElecDensity = val; << 627 } << 628 } 592 } 629 593 630 //....oooOO0OOooo........oooOO0OOooo........oo 594 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 631 595 632 void G4Material::ComputeDensityEffectOnFly(G4b << 596 size_t G4Material::GetNumberOfMaterials() 633 { 597 { 634 if (! IsLocked()) { << 598 return theMaterialTable.size(); 635 if (nullptr == fIonisation) { << 636 fIonisation = new G4IonisParamMat(this); << 637 } << 638 fIonisation->ComputeDensityEffectOnFly(val << 639 } << 640 } 599 } 641 600 642 //....oooOO0OOooo........oooOO0OOooo........oo 601 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 643 602 644 G4MaterialTable* G4Material::GetMaterialTable( << 603 G4Material* 645 { << 604 G4Material::GetMaterial(const G4String& materialName, G4bool warning) 646 struct Holder { << 605 { 647 G4MaterialTable instance; << 606 // search the material by its name 648 ~Holder() { << 607 for (size_t J=0 ; J<theMaterialTable.size() ; ++J) 649 for(auto item : instance) << 608 { 650 delete item; << 609 if (theMaterialTable[J]->GetName() == materialName) 651 } << 610 { return theMaterialTable[J]; } 652 }; << 653 static Holder _holder; << 654 return &_holder.instance; << 655 } << 656 << 657 //....oooOO0OOooo........oooOO0OOooo........oo << 658 << 659 std::size_t G4Material::GetNumberOfMaterials() << 660 << 661 //....oooOO0OOooo........oooOO0OOooo........oo << 662 << 663 G4Material* G4Material::GetMaterial(const G4St << 664 { << 665 // search the material by its name << 666 for (auto const & j : *GetMaterialTable()) { << 667 if (j->GetName() == materialName) { << 668 return j; << 669 } 611 } 670 } << 612 671 << 672 // the material does not exist in the table 613 // the material does not exist in the table 673 if (warn) { << 614 if (warning) { 674 G4cout << "G4Material::GetMaterial() WARNI << 615 G4cout << "G4Material::GetMaterial() WARNING: The material: " 675 << " does not exist in the table. R << 616 << materialName 676 } << 617 << " does not exist in the table. Return NULL pointer." 677 return nullptr; << 618 << G4endl; >> 619 } >> 620 return 0; 678 } 621 } 679 622 680 //....oooOO0OOooo........oooOO0OOooo........oo 623 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 681 624 682 G4Material* G4Material::GetMaterial(G4double z << 625 G4double G4Material::GetZ() const 683 { << 626 { 684 // search the material by its name << 627 if (fNumberOfElements > 1) { 685 for (auto const & mat : *GetMaterialTable()) << 628 G4cout << "G4Material ERROR in GetZ. The material: " << fName 686 if (1 == mat->GetNumberOfElements() && z = << 629 << " is a mixture."; 687 dens == mat->GetDensity()) << 630 G4Exception ("G4Material::GetZ()", "mat036", FatalException, 688 { << 631 "the Atomic number is not well defined." ); 689 return mat; << 632 } 690 } << 633 return (*theElementVector)[0]->GetZ(); 691 } << 692 return nullptr; << 693 } 634 } 694 635 695 //....oooOO0OOooo........oooOO0OOooo........oo 636 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 696 637 697 G4Material* G4Material::GetMaterial(std::size_ << 638 G4double G4Material::GetA() const 698 { << 639 { 699 // search the material by its name << 640 if (fNumberOfElements > 1) { 700 for (auto const & mat : *GetMaterialTable()) << 641 G4cout << "G4Material ERROR in GetA. The material: " << fName 701 if (nComp == mat->GetNumberOfElements() && << 642 << " is a mixture."; 702 return mat; << 643 G4Exception ("G4Material::GetA()", "mat037", FatalException, 703 } << 644 "the Atomic mass is not well defined." ); 704 } << 645 } 705 return nullptr; << 646 return (*theElementVector)[0]->GetA(); 706 } 647 } 707 648 708 //....oooOO0OOooo........oooOO0OOooo........oo 649 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 709 650 710 G4double G4Material::GetZ() const << 651 G4int G4Material::operator==(const G4Material& right) const 711 { 652 { 712 if (fNumberOfElements > 1) { << 653 return (this == (G4Material *) &right); 713 G4ExceptionDescription ed; << 714 ed << "For material " << fName << " ERROR << 715 << " > 1, which is not allowed"; << 716 G4Exception("G4Material::GetZ()", "mat036" << 717 } << 718 return (*theElementVector)[0]->GetZ(); << 719 } 654 } 720 655 721 //....oooOO0OOooo........oooOO0OOooo........oo 656 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 722 657 723 G4double G4Material::GetA() const << 658 G4int G4Material::operator!=(const G4Material& right) const 724 { 659 { 725 if (fNumberOfElements > 1) { << 660 return (this != (G4Material *) &right); 726 G4ExceptionDescription ed; << 727 ed << "For material " << fName << " ERROR << 728 << " > 1, which is not allowed"; << 729 G4Exception("G4Material::GetA()", "mat036" << 730 } << 731 return (*theElementVector)[0]->GetA(); << 732 } 661 } 733 662 734 //....oooOO0OOooo........oooOO0OOooo........oo 663 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 735 664 736 std::ostream& operator<<(std::ostream& flux, c 665 std::ostream& operator<<(std::ostream& flux, const G4Material* material) 737 { 666 { 738 std::ios::fmtflags mode = flux.flags(); 667 std::ios::fmtflags mode = flux.flags(); 739 flux.setf(std::ios::fixed, std::ios::floatfi << 668 flux.setf(std::ios::fixed,std::ios::floatfield); 740 G4long prec = flux.precision(3); 669 G4long prec = flux.precision(3); 741 << 670 742 flux << " Material: " << std::setw(8) << mat << 671 flux 743 << " " << 672 << " Material: " << std::setw(8) << material->fName 744 << " density: " << std::setw(6) << std << 673 << " " << material->fChemicalFormula << " " 745 << G4BestUnit(material->fDensity, "Volu << 674 << " density: " << std::setw(6) << std::setprecision(3) 746 << std::setprecision(3) << G4BestUnit(m << 675 << G4BestUnit(material->fDensity,"Volumic Mass") 747 << " Nucl.Int.Length: " << std::setw(7 << 676 << " RadL: " << std::setw(7) << std::setprecision(3) 748 << G4BestUnit(material->fNuclInterLen, << 677 << G4BestUnit(material->fRadlen,"Length") 749 << std::setw(30) << " Imean: " << std: << 678 << " Nucl.Int.Length: " << std::setw(7) << std::setprecision(3) 750 << G4BestUnit(material->GetIonisation() << 679 << G4BestUnit(material->fNuclInterLen,"Length") 751 << " temperature: " << std::setw(6) << << 680 << "\n" << std::setw(30) 752 << (material->fTemp) / CLHEP::kelvin << << 681 << " Imean: " << std::setw(7) << std::setprecision(3) 753 << " pressure: " << std::setw(6) << st << 682 << G4BestUnit(material->GetIonisation()->GetMeanExcitationEnergy(), 754 << (material->fPressure) / CLHEP::atmos << 683 "Energy"); 755 << "\n"; << 684 756 << 685 if(material->fState == kStateGas) { 757 for (G4int i = 0; i < material->fNumberOfEle << 686 flux 758 flux << "\n ---> " << (*(material->theEl << 687 << " temperature: " << std::setw(6) << std::setprecision(2) 759 << "\n ElmMassFraction: " << << 688 << (material->fTemp)/kelvin << " K" 760 << (material->fMassFractionVector[i]) << 689 << " pressure: " << std::setw(6) << std::setprecision(2) 761 << " ElmAbundance " << std::setw(6) << 690 << (material->fPressure)/atmosphere << " atm"; 762 << 100 * (material->fVecNbOfAtomsPerV << 691 } 763 << " % \n"; << 692 flux << "\n"; 764 } << 693 765 flux.precision(prec); << 694 for (G4int i=0; i<material->fNumberOfElements; i++) { 766 flux.setf(mode, std::ios::floatfield); << 695 flux 767 << 696 << "\n ---> " << (*(material->theElementVector))[i] 768 if (material->IsExtended()) { << 697 << "\n ElmMassFraction: " 769 static_cast<const G4ExtendedMaterial*>(mat << 698 << std::setw(6)<< std::setprecision(2) 770 } << 699 << (material->fMassFractionVector[i])/perCent << " %" 771 << 700 << " ElmAbundance " << std::setw(6)<< std::setprecision(2) >> 701 << 100*(material->VecNbOfAtomsPerVolume[i]) >> 702 /(material->TotNbOfAtomsPerVolume) >> 703 << " % \n"; >> 704 } >> 705 flux.precision(prec); >> 706 flux.setf(mode,std::ios::floatfield); >> 707 772 return flux; 708 return flux; 773 } 709 } 774 710 775 //....oooOO0OOooo........oooOO0OOooo........oo 711 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 776 712 777 std::ostream& operator<<(std::ostream& flux, c 713 std::ostream& operator<<(std::ostream& flux, const G4Material& material) 778 { 714 { 779 flux << &material; << 715 flux << &material; 780 return flux; 716 return flux; 781 } 717 } 782 718 783 //....oooOO0OOooo........oooOO0OOooo........oo 719 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 784 << 720 785 std::ostream& operator<<(std::ostream& flux, c << 721 std::ostream& operator<<(std::ostream& flux, G4MaterialTable MaterialTable) 786 { 722 { 787 // Dump info for all known materials << 723 //Dump info for all known materials 788 flux << "\n***** Table : Nb of materials = " << 724 flux << "\n***** Table : Nb of materials = " << MaterialTable.size() 789 << 725 << " *****\n" << G4endl; 790 for (auto i : MaterialTable) { << 726 791 flux << i << G4endl << G4endl; << 727 for (size_t i=0; i<MaterialTable.size(); ++i) { >> 728 flux << MaterialTable[i] << G4endl << G4endl; 792 } 729 } 793 730 794 return flux; 731 return flux; 795 } << 732 } 796 << 797 //....oooOO0OOooo........oooOO0OOooo........oo << 798 << 799 G4bool G4Material::IsExtended() const { return << 800 733 801 //....oooOO0OOooo........oooOO0OOooo........oo 734 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 802 << 803 void G4Material::SetMaterialPropertiesTable(G4 << 804 { << 805 if (fMaterialPropertiesTable != anMPT && ! I << 806 delete fMaterialPropertiesTable; << 807 fMaterialPropertiesTable = anMPT; << 808 } << 809 } << 810 << 811 //....oooOO0OOooo........oooOO0OOooo........oo << 812 << 813 G4bool G4Material::IsLocked() << 814 { << 815 auto state = G4StateManager::GetStateManager << 816 return state != G4State_PreInit && state != << 817 } << 818 735