source: trunk/source/materials/src/G4Material.cc@ 822

Last change on this file since 822 was 822, checked in by garnier, 17 years ago

import all except CVS
File size: 21.1 KB
Line 
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//
27// $Id: G4Material.cc,v 1.38.2.1 2008/04/28 07:24:36 gcosmo Exp $
28// GEANT4 tag $Name: geant4-09-01-patch-02 $
29//
30//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
31//
32// 26-06-96, Code uses operators (+=, *=, ++, -> etc.) correctly, P. Urban
33// 10-07-96, new data members added by L.Urban
34// 12-12-96, new data members added by L.Urban
35// 20-01-97, aesthetic rearrangement. RadLength calculation modified.
36// Data members Zeff and Aeff REMOVED (i.e. passed to the Elements).
37// (local definition of Zeff in DensityEffect and FluctModel...)
38// Vacuum defined as a G4State. Mixture flag removed, M.Maire.
39// 29-01-97, State=Vacuum automatically set density=0 in the contructors.
40// Subsequent protections have been put in the calculation of
41// MeanExcEnergy, ShellCorrectionVector, DensityEffect, M.Maire.
42// 11-02-97, ComputeDensityEffect() rearranged, M.Maire.
43// 20-03-97, corrected initialization of pointers, M.Maire.
44// 28-05-98, the kState=kVacuum has been removed.
45// automatic check for a minimal density, M.Maire
46// 12-06-98, new method AddMaterial() allowing mixture of materials, M.Maire
47// 09-07-98, ionisation parameters removed from the class, M.Maire
48// 05-10-98, change names: NumDensity -> NbOfAtomsPerVolume
49// 18-11-98, new interface to SandiaTable
50// 19-01-99 enlarge tolerance on test of coherence of gas conditions
51// 19-07-99, Constructors with chemicalFormula added by V.Ivanchenko
52// 16-01-01, Nuclear interaction length, M.Maire
53// 12-03-01, G4bool fImplicitElement;
54// copy constructor and assignement operator revised (mma)
55// 03-05-01, flux.precision(prec) at begin/end of operator<<
56// 17-07-01, migration to STL. M. Verderi.
57// 14-09-01, Suppression of the data member fIndexInTable
58// 26-02-02, fIndexInTable renewed
59// 16-04-02, G4Exception put in constructor with chemical formula
60// 06-05-02, remove the check of the ideal gas state equation
61// 06-08-02, remove constructors with chemical formula (mma)
62// 22-01-04, proper STL handling of theElementVector (Hisaya)
63// 30-03-05, warning in GetMaterial(materialName)
64// 09-03-06, minor change of printout (V.Ivanchenko)
65// 10-01-07, compute fAtomVector in the case of mass fraction (V.Ivanchenko)
66// 27-07-07, improve destructor (V.Ivanchenko)
67// 18-10-07, move definition of material index to InitialisePointers (V.Ivanchenko)
68//
69
70//
71//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
72
73#include "G4Material.hh"
74#include "G4UnitsTable.hh"
75#include <iomanip>
76
77
78G4MaterialTable G4Material::theMaterialTable;
79
80//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
81
82// Constructor to create a material from scratch
83
84G4Material::G4Material(const G4String& name, G4double z,
85 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 { G4cerr << "--- Warning from G4Material::G4Material()"
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)
96 << "g/cm3" << G4endl;
97 density = universe_mean_density;
98 }
99
100 fDensity = density;
101 fState = state;
102 fTemp = temp;
103 fPressure = pressure;
104 fChemicalFormula = " ";
105
106 // Initialize theElementVector allocating one
107 // element corresponding to this material
108 maxNbComponents = fNumberOfComponents = fNumberOfElements = 1;
109 fImplicitElement = true;
110 theElementVector = new G4ElementVector();
111 theElementVector->push_back( new G4Element(name, " ", z, a));
112 fMassFractionVector = new G4double[1];
113 fMassFractionVector[0] = 1. ;
114
115 (*theElementVector)[0] -> increaseCountUse();
116
117 if (fState == kStateUndefined)
118 {
119 if (fDensity > kGasThreshold) fState = kStateSolid;
120 else fState = kStateGas;
121 }
122
123 ComputeDerivedQuantities();
124}
125
126//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
127
128// Constructor to create a material from a List of constituents
129// (elements and/or materials) added with AddElement or AddMaterial
130
131G4Material::G4Material(const G4String& name, G4double density,
132 G4int nComponents,
133 G4State state, G4double temp, G4double pressure)
134 : fName(name)
135{
136 InitializePointers();
137
138 if (density < universe_mean_density)
139 {G4cerr << "--- Warning from G4Material::G4Material()"
140 << " define a material with density=0 is not allowed. \n"
141 << " The material " << name << " will be constructed with the"
142 << " default minimal density: " << universe_mean_density/(g/cm3)
143 << "g/cm3" << G4endl;
144 density = universe_mean_density;
145 }
146
147 fDensity = density;
148 fState = state;
149 fTemp = temp;
150 fPressure = pressure;
151 fChemicalFormula = " ";
152
153 maxNbComponents = nComponents;
154 fNumberOfComponents = fNumberOfElements = 0;
155 fImplicitElement = false;
156 theElementVector = new G4ElementVector();
157 theElementVector->reserve(maxNbComponents);
158
159 if (fState == kStateUndefined)
160 {
161 if (fDensity > kGasThreshold) fState = kStateSolid;
162 else fState = kStateGas;
163 }
164}
165
166//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
167
168// Fake default constructor - sets only member data and allocates memory
169// for usage restricted to object persistency
170
171G4Material::G4Material(__void__&)
172 : fNumberOfComponents(0), fNumberOfElements(0), theElementVector(0),
173 fImplicitElement(false), fMassFractionVector(0), fAtomsVector(0),
174 fMaterialPropertiesTable(0), fIndexInTable(0),
175 VecNbOfAtomsPerVolume(0), fIonisation(0), fSandiaTable(0)
176{
177}
178
179//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
180
181// AddElement -- composition by atom count
182
183void G4Material::AddElement(G4Element* element, G4int nAtoms)
184{
185 // initialization
186 if ( fNumberOfElements == 0 ) {
187 fAtomsVector = new G4int [maxNbComponents];
188 fMassFractionVector = new G4double[maxNbComponents];
189 }
190
191 // filling ...
192 if ( G4int(fNumberOfElements) < maxNbComponents ) {
193 theElementVector->push_back(element);
194 fAtomsVector [fNumberOfElements] = nAtoms;
195 fNumberOfComponents = ++fNumberOfElements;
196 element->increaseCountUse();
197 }
198 else
199 G4Exception
200 ("ERROR!!! - Attempt to add more than the declared number of elements.");
201
202 // filled.
203 if ( G4int(fNumberOfElements) == maxNbComponents ) {
204 // compute proportion by mass
205 size_t i=0;
206 G4double Zmol(0.), Amol(0.);
207 for (i=0;i<fNumberOfElements;i++) {
208 Zmol += fAtomsVector[i]*(*theElementVector)[i]->GetZ();
209 Amol += fAtomsVector[i]*(*theElementVector)[i]->GetA();
210 }
211 for (i=0;i<fNumberOfElements;i++) {
212 fMassFractionVector[i] = fAtomsVector[i]
213 *(*theElementVector)[i]->GetA()/Amol;
214 }
215
216 ComputeDerivedQuantities();
217 }
218}
219
220//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
221
222// AddElement -- composition by fraction of mass
223
224void G4Material::AddElement(G4Element* element, G4double fraction)
225{
226 // initialization
227 if (fNumberOfComponents == 0) {
228 fMassFractionVector = new G4double[50];
229 fAtomsVector = new G4int [50];
230 }
231
232 // filling ...
233 if (G4int(fNumberOfComponents) < maxNbComponents) {
234 size_t el = 0;
235 while ((el<fNumberOfElements)&&(element!=(*theElementVector)[el])) el++;
236 if (el<fNumberOfElements) fMassFractionVector[el] += fraction;
237 else {
238 theElementVector->push_back(element);
239 fMassFractionVector[el] = fraction;
240 fNumberOfElements++;
241 element->increaseCountUse();
242 }
243 fNumberOfComponents++;
244 }
245 else
246 G4Exception
247 ("ERROR!!! - Attempt to add more than the declared number of components.");
248
249 // filled.
250 if (G4int(fNumberOfComponents) == maxNbComponents) {
251
252 size_t i=0;
253 G4double Zmol(0.), Amol(0.);
254 // check sum of weights -- OK?
255 G4double wtSum(0.0);
256 for (i=0;i<fNumberOfElements;i++) {
257 wtSum += fMassFractionVector[i];
258 Zmol += fMassFractionVector[i]*(*theElementVector)[i]->GetZ();
259 Amol += fMassFractionVector[i]*(*theElementVector)[i]->GetA();
260 }
261 if (std::abs(1.-wtSum) > perThousand) {
262 G4cerr << "WARNING !! for " << fName << " sum of fractional masses "
263 << wtSum << " is not 1 - results may be wrong"
264 << G4endl;
265 }
266 for (i=0;i<fNumberOfElements;i++) {
267 fAtomsVector[i] =
268 G4int(fMassFractionVector[i]*Amol/(*theElementVector)[i]->GetA()+0.5);
269 }
270
271 ComputeDerivedQuantities();
272 }
273}
274
275//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
276
277// AddMaterial -- composition by fraction of mass
278
279void G4Material::AddMaterial(G4Material* material, G4double fraction)
280{
281 // initialization
282 if (fNumberOfComponents == 0) {
283 fMassFractionVector = new G4double[50];
284 fAtomsVector = new G4int [50];
285 }
286
287 // filling ...
288 if (G4int(fNumberOfComponents) < maxNbComponents) {
289 for (size_t elm=0; elm < material->GetNumberOfElements(); elm++)
290 {
291 G4Element* element = (*(material->GetElementVector()))[elm];
292 size_t el = 0;
293 while ((el<fNumberOfElements)&&(element!=(*theElementVector)[el])) el++;
294 if (el < fNumberOfElements) fMassFractionVector[el] += fraction
295 *(material->GetFractionVector())[elm];
296 else {
297 theElementVector->push_back(element);
298 fMassFractionVector[el] = fraction
299 *(material->GetFractionVector())[elm];
300 fNumberOfElements++;
301 element->increaseCountUse();
302 }
303 }
304 fNumberOfComponents++;
305 }
306 else
307 G4Exception
308 ("ERROR!!! - Attempt to add more than the declared number of components.");
309
310 // filled.
311 if (G4int(fNumberOfComponents) == maxNbComponents) {
312 size_t i=0;
313 G4double Zmol(0.), Amol(0.);
314 // check sum of weights -- OK?
315 G4double wtSum(0.0);
316 for (i=0;i<fNumberOfElements;i++) {
317 wtSum += fMassFractionVector[i];
318 Zmol += fMassFractionVector[i]*(*theElementVector)[i]->GetZ();
319 Amol += fMassFractionVector[i]*(*theElementVector)[i]->GetA();
320 }
321 if (std::abs(1.-wtSum) > perThousand) {
322 G4cerr << "WARNING !! for " << fName << " sum of fractional masses "
323 << wtSum << " is not 1 - results may be wrong"
324 << G4endl;
325 }
326 for (i=0;i<fNumberOfElements;i++) {
327 fAtomsVector[i] =
328 G4int(fMassFractionVector[i]*Amol/(*theElementVector)[i]->GetA()+0.5);
329 }
330
331 ComputeDerivedQuantities();
332 }
333}
334
335//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
336
337void G4Material::ComputeDerivedQuantities()
338{
339 // Header routine to compute various properties of material.
340 //
341
342 // Number of atoms per volume (per element), total nb of electrons per volume
343 G4double Zi, Ai;
344 TotNbOfAtomsPerVolume = 0.;
345 if (VecNbOfAtomsPerVolume) delete [] VecNbOfAtomsPerVolume;
346 VecNbOfAtomsPerVolume = new G4double[fNumberOfElements];
347 TotNbOfElectPerVolume = 0.;
348 for (size_t i=0;i<fNumberOfElements;i++) {
349 Zi = (*theElementVector)[i]->GetZ();
350 Ai = (*theElementVector)[i]->GetA();
351 VecNbOfAtomsPerVolume[i] = Avogadro*fDensity*fMassFractionVector[i]/Ai;
352 TotNbOfAtomsPerVolume += VecNbOfAtomsPerVolume[i];
353 TotNbOfElectPerVolume += VecNbOfAtomsPerVolume[i]*Zi;
354 }
355
356 ComputeRadiationLength();
357 ComputeNuclearInterLength();
358
359 if (fIonisation) delete fIonisation;
360 fIonisation = new G4IonisParamMat(this);
361 if (fSandiaTable) delete fSandiaTable;
362 fSandiaTable = new G4SandiaTable(this);
363}
364
365//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
366
367void G4Material::ComputeRadiationLength()
368{
369 G4double radinv = 0.0 ;
370 for (size_t i=0;i<fNumberOfElements;i++) {
371 radinv += VecNbOfAtomsPerVolume[i]*((*theElementVector)[i]->GetfRadTsai());
372 }
373 fRadlen = (radinv <= 0.0 ? DBL_MAX : 1./radinv);
374}
375
376//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
377
378void G4Material::ComputeNuclearInterLength()
379{
380 const G4double lambda0 = 35*g/cm2;
381 G4double NILinv = 0.0;
382 for (size_t i=0;i<fNumberOfElements;i++) {
383 NILinv +=
384 VecNbOfAtomsPerVolume[i]*std::pow(((*theElementVector)[i]->GetN()),0.6666667);
385 }
386 NILinv *= amu/lambda0;
387 fNuclInterLen = (NILinv <= 0.0 ? DBL_MAX : 1./NILinv);
388}
389
390//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
391
392void G4Material::InitializePointers()
393{
394 theElementVector = 0;
395 fMassFractionVector = 0;
396 fAtomsVector = 0;
397 fMaterialPropertiesTable = 0;
398
399 VecNbOfAtomsPerVolume = 0;
400 fIonisation = 0;
401 fSandiaTable = 0;
402
403 // Store in the static Table of Materials
404 theMaterialTable.push_back(this);
405 fIndexInTable = theMaterialTable.size() - 1;
406}
407
408//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
409
410const G4MaterialTable* G4Material::GetMaterialTable()
411{
412 return &theMaterialTable;
413}
414
415//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
416
417size_t G4Material::GetNumberOfMaterials()
418{
419 return theMaterialTable.size();
420}
421
422//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
423
424G4Material* G4Material::GetMaterial(G4String materialName, G4bool warning)
425{
426 // search the material by its name
427 for (size_t J=0 ; J<theMaterialTable.size() ; J++)
428 {
429 if (theMaterialTable[J]->GetName() == materialName)
430 return theMaterialTable[J];
431 }
432
433 // the material does not exist in the table
434 if (warning) {
435 G4cout << "\n---> warning from G4Material::GetMaterial(). The material: "
436 << materialName << " does not exist in the table. Return NULL pointer."
437 << G4endl;
438 }
439 return 0;
440}
441
442//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
443
444G4Material::G4Material(const G4Material& right)
445{
446 InitializePointers();
447 *this = right;
448}
449
450//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
451
452G4Material::~G4Material()
453{
454 // G4cout << "### Destruction of material " << fName << " started" <<G4endl;
455 if (theElementVector) delete theElementVector;
456 if (fMassFractionVector) delete [] fMassFractionVector;
457 if (fAtomsVector) delete [] fAtomsVector;
458 if (VecNbOfAtomsPerVolume) delete [] VecNbOfAtomsPerVolume;
459 if (fIonisation) delete fIonisation;
460 if (fSandiaTable) delete fSandiaTable;
461
462 // Remove this material from theMaterialTable.
463 //
464 theMaterialTable[fIndexInTable] = 0;
465}
466
467//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
468
469const G4Material& G4Material::operator=(const G4Material& right)
470{
471 if (this != &right)
472 {
473 fName = right.fName;
474 fChemicalFormula = right.fChemicalFormula;
475 fDensity = right.fDensity;
476 fState = right.fState;
477 fTemp = right.fTemp;
478 fPressure = right.fPressure;
479
480 if (fImplicitElement) delete ((*theElementVector)[0]);
481 if (theElementVector) delete theElementVector;
482 if (fMassFractionVector) delete [] fMassFractionVector;
483 if (fAtomsVector) delete [] fAtomsVector;
484
485 maxNbComponents = right.maxNbComponents;
486 fNumberOfComponents = right.fNumberOfComponents;
487 fNumberOfElements = right.fNumberOfElements;
488 fImplicitElement = right.fImplicitElement;
489
490 if (fImplicitElement) {
491 G4double z = (*right.theElementVector)[0]->GetZ();
492 G4double a = (*right.theElementVector)[0]->GetA();
493 theElementVector = new G4ElementVector(1,(G4Element*)0);
494 (*theElementVector)[0] = new G4Element(fName," ",z,a);
495 fMassFractionVector = new G4double[1];
496 fMassFractionVector[0] = 1.;
497 } else {
498 theElementVector = new G4ElementVector(fNumberOfElements,0);
499 fMassFractionVector = new G4double[fNumberOfElements];
500 for (size_t i=0; i<fNumberOfElements; i++) {
501 (*theElementVector)[i]= (*right.theElementVector)[i];
502 fMassFractionVector[i]= right.fMassFractionVector[i];
503 }
504 }
505
506 if (right.fAtomsVector) {
507 fAtomsVector = new G4int[fNumberOfElements];
508 for (size_t i=0; i<fNumberOfElements; i++)
509 fAtomsVector[i] = right.fAtomsVector[i];
510 }
511
512 fMaterialPropertiesTable = right.fMaterialPropertiesTable;
513
514 ComputeDerivedQuantities();
515 }
516 return *this;
517}
518
519//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
520
521G4int G4Material::operator==(const G4Material& right) const
522{
523 return (this == (G4Material *) &right);
524}
525
526//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
527
528G4int G4Material::operator!=(const G4Material& right) const
529{
530 return (this != (G4Material *) &right);
531}
532
533//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
534
535
536std::ostream& operator<<(std::ostream& flux, G4Material* material)
537{
538 std::ios::fmtflags mode = flux.flags();
539 flux.setf(std::ios::fixed,std::ios::floatfield);
540 G4long prec = flux.precision(3);
541
542 flux
543 << " Material: " << std::setw(8) << material->fName
544 << " " << material->fChemicalFormula << " "
545 << " density: " << std::setw(6) << std::setprecision(3)
546 << G4BestUnit(material->fDensity,"Volumic Mass")
547 << " RadL: " << std::setw(7) << std::setprecision(3)
548 << G4BestUnit(material->fRadlen,"Length")
549 << " Imean: " << std::setw(7) << std::setprecision(3)
550 << G4BestUnit(material->GetIonisation()->GetMeanExcitationEnergy(),"Energy");
551 if(material->fState == kStateGas)
552 flux
553 << " temperature: " << std::setw(6) << std::setprecision(2)
554 << (material->fTemp)/kelvin << " K"
555 << " pressure: " << std::setw(6) << std::setprecision(2)
556 << (material->fPressure)/atmosphere << " atm";
557
558 for (size_t i=0; i<material->fNumberOfElements; i++)
559 flux
560 << "\n ---> " << (*(material->theElementVector))[i]
561 << " ElmMassFraction: " << std::setw(6)<< std::setprecision(2)
562 << (material->fMassFractionVector[i])/perCent << " %"
563 << " ElmAbundance " << std::setw(6)<< std::setprecision(2)
564 << 100*(material->VecNbOfAtomsPerVolume[i])/(material->TotNbOfAtomsPerVolume)
565 << " %";
566
567 flux.precision(prec);
568 flux.setf(mode,std::ios::floatfield);
569
570 return flux;
571}
572
573//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
574
575 std::ostream& operator<<(std::ostream& flux, G4Material& material)
576{
577 flux << &material;
578 return flux;
579}
580
581//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
582
583std::ostream& operator<<(std::ostream& flux, G4MaterialTable MaterialTable)
584{
585 //Dump info for all known materials
586 flux << "\n***** Table : Nb of materials = " << MaterialTable.size()
587 << " *****\n" << G4endl;
588
589 for (size_t i=0; i<MaterialTable.size(); i++) flux << MaterialTable[i]
590 << G4endl << G4endl;
591
592 return flux;
593}
594
595//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
Note: See TracBrowser for help on using the repository browser.