Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/materials/src/G4Material.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

Diff markup

Differences between /materials/src/G4Material.cc (Version 11.3.0) and /materials/src/G4Material.cc (Version 10.0)


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