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