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

Last change on this file since 843 was 822, checked in by garnier, 16 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.