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 5.2.p2)


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