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