Index: trunk/source/materials/src/G4AtomicShells.cc
===================================================================
--- trunk/source/materials/src/G4AtomicShells.cc	(revision 822)
+++ trunk/source/materials/src/G4AtomicShells.cc	(revision 850)
@@ -26,5 +26,5 @@
 //
 // $Id: G4AtomicShells.cc,v 1.7 2006/10/17 15:15:46 vnivanch Exp $
-// GEANT4 tag $Name:  $
+// GEANT4 tag $Name: HEAD $
 //
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
Index: trunk/source/materials/src/G4Element.cc
===================================================================
--- trunk/source/materials/src/G4Element.cc	(revision 822)
+++ trunk/source/materials/src/G4Element.cc	(revision 850)
@@ -25,6 +25,6 @@
 //
 //
-// $Id: G4Element.cc,v 1.26 2007/10/18 11:14:33 vnivanch Exp $
-// GEANT4 tag $Name:  $
+// $Id: G4Element.cc,v 1.31 2008/08/11 11:53:11 vnivanch Exp $
+// GEANT4 tag $Name: HEAD $
 //
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
@@ -69,14 +69,18 @@
   : fName(name), fSymbol(symbol)		     
 {
-  if (zeff<1.) G4Exception (" ERROR from G4Element::G4Element !"
-			    " It is not allowed to create an Element with Z < 1" );
-
-  if (aeff/(g/mole)<zeff) G4Exception (" ERROR from G4Element::G4Element !"
-				       " Attempt to create an Element with N < Z !!!" );
-
-  if ((zeff-G4int(zeff)) > perMillion)
+  G4int iz = (G4int)zeff;
+  if (zeff<1.) {
+    G4cout << "G4Element ERROR:  " << name << " Z= " << zeff 
+	   << " A= " << aeff/(g/mole) << G4endl; 
+    G4Exception (" ERROR from G4Element::G4Element !"
+		 " It is not allowed to create an Element with Z < 1" );
+  }
+  if (std::abs(zeff - iz) > perMillion) {
+    G4cout << "G4Element Warning:  " << name << " Z= " << zeff 
+	   << " A= " << aeff/(g/mole) << G4endl; 
     G4cerr << name << " : WARNING from G4Element::G4Element !"  
       " Trying to define an element as a mixture directly via effective Z."
 	   << G4endl;
+  }
 
   InitializePointers();
@@ -85,10 +89,19 @@
   fNeff   = aeff/(g/mole);
   fAeff   = aeff;
+
+  if(fNeff < 1.0) fNeff = 1.0;
+
+  if (fNeff < zeff) {
+    G4cout << "G4Element ERROR:  " << name << " Z= " << zeff 
+	   << " A= " << fNeff << G4endl; 
+    G4Exception (" ERROR from G4Element::G4Element !"
+		 " Attempt to create an Element with N < Z !!!" );
+  }
    
-  fNbOfAtomicShells = G4AtomicShells::GetNumberOfShells((G4int)fZeff);
+  fNbOfAtomicShells = G4AtomicShells::GetNumberOfShells(iz);
   fAtomicShells     = new G4double[fNbOfAtomicShells];
-  for (G4int i=0;i<fNbOfAtomicShells;i++)
-    fAtomicShells[i] = G4AtomicShells::GetBindingEnergy((G4int)fZeff,i);
-           
+  for (G4int i=0;i<fNbOfAtomicShells;i++) {
+    fAtomicShells[i] = G4AtomicShells::GetBindingEnergy(iz, i);
+  }
   ComputeDerivedQuantities();
 }
@@ -117,15 +130,19 @@
 void G4Element::AddIsotope(G4Isotope* isotope, G4double abundance)
 {
-  if (theIsotopeVector == 0)
+  if (theIsotopeVector == 0) {
+    G4cout << "G4Element ERROR:  " << fName << G4endl; 
     G4Exception ("ERROR from G4Element::AddIsotope!"
 		 " Trying to add an Isotope before contructing the element.");
+  }
+  G4int iz = isotope->GetZ();
 
   // filling ...
   if ( fNumberOfIsotopes < theIsotopeVector->size() ) {
     // check same Z
-    if (fNumberOfIsotopes==0) fZeff = G4double(isotope->GetZ());
-    else if (G4double(isotope->GetZ()) != fZeff) 
+    if (fNumberOfIsotopes==0) fZeff = G4double(iz);
+    else if (G4double(iz) != fZeff) { 
       G4Exception ("ERROR from G4Element::AddIsotope!"
 		   " Try to add isotopes with different Z");
+    }
     //Z ok   
     fRelativeAbundanceVector[fNumberOfIsotopes] = abundance;
@@ -133,7 +150,9 @@
     ++fNumberOfIsotopes;
     isotope->increaseCountUse();
-  } 
-  else G4Exception ("ERROR from G4Element::AddIsotope!"  
-		    " Attempt to add more than the declared number of isotopes.");
+  } else {
+    G4cout << "G4Element ERROR:  " << fName << G4endl; 
+    G4Exception ("ERROR from G4Element::AddIsotope!"  
+		 " Attempt to add more than the declared number of isotopes.");
+  }
 
   // filled.
@@ -151,8 +170,9 @@
     fAeff /=  wtSum;
       
-    fNbOfAtomicShells = G4AtomicShells::GetNumberOfShells((G4int)fZeff);
+    fNbOfAtomicShells = G4AtomicShells::GetNumberOfShells(iz);
     fAtomicShells     = new G4double[fNbOfAtomicShells];
-    for (G4int j=0;j<fNbOfAtomicShells;j++)
-      fAtomicShells[j] = G4AtomicShells::GetBindingEnergy((G4int)fZeff,j);
+    for (G4int j=0;j<fNbOfAtomicShells;j++) {
+      fAtomicShells[j] = G4AtomicShells::GetBindingEnergy(iz, j);
+    }
          
     ComputeDerivedQuantities();
@@ -213,7 +233,7 @@
   // check if elements with same Z already exist
   fIndexZ = 0; 
-  for (size_t J=0 ; J<fIndexInTable ; J++)
+  for (size_t J=0 ; J<fIndexInTable ; J++) {
     if (theElementTable[J]->GetZ() == fZeff) fIndexZ++;
-
+  }
   //nb of materials which use this element
   fCountUse = 0;
@@ -257,5 +277,5 @@
 
   G4double Lrad, Lprad;
-  G4int iz = (int)(fZeff+0.5) - 1 ;
+  G4int iz = (G4int)(fZeff+0.5) - 1 ;
   if (iz <= 3) { Lrad = Lrad_light[iz] ;  Lprad = Lprad_light[iz] ; }
     else { Lrad = std::log(184.15) - logZ3 ; Lprad = std::log(1194.) - 2*logZ3;}
Index: trunk/source/materials/src/G4IonisParamElm.cc
===================================================================
--- trunk/source/materials/src/G4IonisParamElm.cc	(revision 822)
+++ trunk/source/materials/src/G4IonisParamElm.cc	(revision 850)
@@ -25,6 +25,6 @@
 //
 //
-// $Id: G4IonisParamElm.cc,v 1.14 2006/06/29 19:12:43 gunter Exp $
-// GEANT4 tag $Name:  $
+// $Id: G4IonisParamElm.cc,v 1.15 2008/06/03 14:30:10 vnivanch Exp $
+// GEANT4 tag $Name: HEAD $
 //
 //
@@ -47,19 +47,19 @@
 G4IonisParamElm::G4IonisParamElm(G4double Z)
 {
-    if (Z < 1.) G4Exception
-      (" ERROR! It is not allowed to create an Element with Z < 1" );
+  if (Z < 1.) 
+    G4Exception (" ERROR! It is not allowed to create an Element with Z < 1" );
     
-    // some basic functions of the atomic number
-    fZ = Z;
-    fZ3  = std::pow(fZ, 1./3.);
-    fZZ3 = std::pow(fZ*(fZ+1.), 1./3.);
-    flogZ3 = std::log(fZ)/3.;
+  // some basic functions of the atomic number
+  fZ = Z;
+  fZ3  = std::pow(fZ, 1./3.);
+  fZZ3 = std::pow(fZ*(fZ+1.), 1./3.);
+  flogZ3 = std::log(fZ)/3.;
    
-     // Parameters for energy loss by ionisation   
-
-     // Mean excitation energy
-     // from "Stopping Powers for Electrons and Positrons"
-     // ICRU Report N#37, 1984  (energy in eV)
-    static double exc[100] = { 
+  // Parameters for energy loss by ionisation   
+
+  // Mean excitation energy
+  // from "Stopping Powers for Electrons and Positrons"
+  // ICRU Report N#37, 1984  (energy in eV)
+  static double exc[100] = { 
     21.8, 20.9,	13.3, 15.9, 15.2, 13.0,	11.7, 11.9, 11.5, 13.7,
     13.6, 13.0,	12.8, 12.4, 11.5, 11.3,	10.2, 10.4, 10.0,  9.6,
@@ -73,28 +73,64 @@
      9.6,  9.7,  9.7,  9.8,  9.8,  9.8,  9.8,  9.9,  9.9,  9.9 };
 
-    G4int iz = (G4int)Z - 1 ;
-    if(0  > iz) iz = 0;
-    else if(99 < iz) iz = 99 ;
-    fMeanExcitationEnergy = fZ * exc[iz] * eV ;
-
-
-    fTau0 = 0.1*fZ3*MeV/proton_mass_c2;
-    fTaul = 2.*MeV/proton_mass_c2;
-
-    // compute the Bethe-Bloch formula for energy = fTaul*particle mass
-    G4double rate = fMeanExcitationEnergy/electron_mass_c2 ;
-    G4double w = fTaul*(fTaul+2.) ;
-    fBetheBlochLow = (fTaul+1.)*(fTaul+1.)*std::log(2.*w/rate)/w - 1. ;
-    fBetheBlochLow = 2.*fZ*twopi_mc2_rcl2*fBetheBlochLow ; 
+  G4int iz = (G4int)Z - 1 ;
+  if(0  > iz) iz = 0;
+  else if(99 < iz) iz = 99 ;
+  fMeanExcitationEnergy = fZ * exc[iz] * eV ;
+
+  // compute parameters for ion transport
+  // The aproximation from:
+  // J.F.Ziegler, J.P. Biersack, U. Littmark
+  // The Stopping and Range of Ions in Matter,
+  // Vol.1, Pergamon Press, 1985
+  // Fast ions or hadrons
+
+  if(91 < iz) iz = 91;
+
+  static G4double vFermi[92] = {
+    1.0309,  0.15976, 0.59782, 1.0781,  1.0486,  1.0,     1.058,   0.93942, 0.74562, 0.3424,
+    0.45259, 0.71074, 0.90519, 0.97411, 0.97184, 0.89852, 0.70827, 0.39816, 0.36552, 0.62712,
+    0.81707, 0.9943,  1.1423,  1.2381,  1.1222,  0.92705, 1.0047,  1.2,     1.0661,  0.97411,
+    0.84912, 0.95,    1.0903,  1.0429,  0.49715, 0.37755, 0.35211, 0.57801, 0.77773, 1.0207,
+    1.029,   1.2542,  1.122,   1.1241,  1.0882,  1.2709,  1.2542,  0.90094, 0.74093, 0.86054,
+    0.93155, 1.0047,  0.55379, 0.43289, 0.32636, 0.5131,  0.695,   0.72591, 0.71202, 0.67413,
+    0.71418, 0.71453, 0.5911,  0.70263, 0.68049, 0.68203, 0.68121, 0.68532, 0.68715, 0.61884,
+    0.71801, 0.83048, 1.1222,  1.2381,  1.045,   1.0733,  1.0953,  1.2381,  1.2879,  0.78654,
+    0.66401, 0.84912, 0.88433, 0.80746, 0.43357, 0.41923, 0.43638, 0.51464, 0.73087, 0.81065,
+    1.9578,  1.0257} ;
+
+  static G4double lFactor[92] = {
+    1.0,  1.0,  1.1,  1.06, 1.01, 1.03, 1.04, 0.99, 0.95, 0.9,
+    0.82, 0.81, 0.83, 0.88, 1.0,  0.95, 0.97, 0.99, 0.98, 0.97,
+    0.98, 0.97, 0.96, 0.93, 0.91, 0.9,  0.88, 0.9,  0.9,  0.9,
+    0.9,  0.85, 0.9,  0.9,  0.91, 0.92, 0.9,  0.9,  0.9,  0.9,
+    0.9,  0.88, 0.9,  0.88, 0.88, 0.9,  0.9,  0.88, 0.9,  0.9,
+    0.9,  0.9,  0.96, 1.2,  0.9,  0.88, 0.88, 0.85, 0.9,  0.9,
+    0.92, 0.95, 0.99, 1.03, 1.05, 1.07, 1.08, 1.1,  1.08, 1.08,
+    1.08, 1.08, 1.09, 1.09, 1.1,  1.11, 1.12, 1.13, 1.14, 1.15,
+    1.17, 1.2,  1.18, 1.17, 1.17, 1.16, 1.16, 1.16, 1.16, 1.16,
+    1.16, 1.16} ;
+
+  fVFermi  = vFermi[iz];
+  fLFactor = lFactor[iz];
+
+  // obsolete parameters for ionisation
+  fTau0 = 0.1*fZ3*MeV/proton_mass_c2;
+  fTaul = 2.*MeV/proton_mass_c2;
+
+  // compute the Bethe-Bloch formula for energy = fTaul*particle mass
+  G4double rate = fMeanExcitationEnergy/electron_mass_c2 ;
+  G4double w = fTaul*(fTaul+2.) ;
+  fBetheBlochLow = (fTaul+1.)*(fTaul+1.)*std::log(2.*w/rate)/w - 1. ;
+  fBetheBlochLow = 2.*fZ*twopi_mc2_rcl2*fBetheBlochLow ; 
   
-    fClow = std::sqrt(fTaul)*fBetheBlochLow;
-    fAlow = 6.458040 * fClow/fTau0;
-    G4double Taum = 0.035*fZ3*MeV/proton_mass_c2;
-    fBlow =-3.229020*fClow/(fTau0*std::sqrt(Taum));
-
-    // Shell correction factors
-    fShellCorrectionVector = new G4double[3]; //[3]
-    rate = 0.001*fMeanExcitationEnergy/eV;
-    G4double rate2 = rate*rate;
+  fClow = std::sqrt(fTaul)*fBetheBlochLow;
+  fAlow = 6.458040 * fClow/fTau0;
+  G4double Taum = 0.035*fZ3*MeV/proton_mass_c2;
+  fBlow =-3.229020*fClow/(fTau0*std::sqrt(Taum));
+
+  // Shell correction parameterization
+  fShellCorrectionVector = new G4double[3]; //[3]
+  rate = 0.001*fMeanExcitationEnergy/eV;
+  G4double rate2 = rate*rate;
     /*    
     fShellCorrectionVector[0] = ( 1.10289e5 + 5.14781e8*rate)*rate2 ;
@@ -102,7 +138,7 @@
     fShellCorrectionVector[2] = (-9.92256e1 + 2.10823e5*rate)*rate2 ;
     */
-    fShellCorrectionVector[0] = ( 0.422377   + 3.858019*rate)*rate2 ;
-    fShellCorrectionVector[1] = ( 0.0304043  - 0.1667989*rate)*rate2 ;
-    fShellCorrectionVector[2] = (-0.00038106 + 0.00157955*rate)*rate2 ;
+  fShellCorrectionVector[0] = ( 0.422377   + 3.858019*rate)*rate2 ;
+  fShellCorrectionVector[1] = ( 0.0304043  - 0.1667989*rate)*rate2 ;
+  fShellCorrectionVector[2] = (-0.00038106 + 0.00157955*rate)*rate2 ;
 }
 
Index: trunk/source/materials/src/G4IonisParamMat.cc
===================================================================
--- trunk/source/materials/src/G4IonisParamMat.cc	(revision 822)
+++ trunk/source/materials/src/G4IonisParamMat.cc	(revision 850)
@@ -25,6 +25,6 @@
 //
 //
-// $Id: G4IonisParamMat.cc,v 1.20 2007/09/27 14:05:47 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-01-patch-02 $
+// $Id: G4IonisParamMat.cc,v 1.25 2008/07/08 10:34:56 vnivanch Exp $
+// GEANT4 tag $Name: HEAD $
 //
 // 
@@ -39,4 +39,5 @@
 // 10-05-05, add a missing coma in FindMeanExcitationEnergy() - Bug#746 (mma)
 // 27-09-07, add computation of parameters for ions (V.Ivanchenko)
+// 04-03-08, remove reference to G4NistManager. Add fBirks constant (mma)
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
@@ -44,5 +45,4 @@
 #include "G4IonisParamMat.hh"
 #include "G4Material.hh"
-#include "G4NistManager.hh"
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
@@ -55,4 +55,7 @@
   ComputeFluctModel();
   ComputeIonParameters();
+  
+  fBirks = 0.;
+  fMeanEnergyPerIon = 0.;
 }
 
@@ -78,15 +81,26 @@
   fLogMeanExcEnergy = 0.;
 
-
-  for (size_t i=0; i < fMaterial->GetNumberOfElements(); i++) {
-    fLogMeanExcEnergy += 
-             (fMaterial->GetVecNbOfAtomsPerVolume())[i]
-            *((*(fMaterial->GetElementVector()))[i]->GetZ())
-            *std::log((*(fMaterial->GetElementVector()))[i]->GetIonisation()
-             ->GetMeanExcitationEnergy());
+  size_t nElements = fMaterial->GetNumberOfElements();
+  const G4ElementVector* elmVector = fMaterial->GetElementVector();
+  const G4double* nAtomsPerVolume = fMaterial->GetVecNbOfAtomsPerVolume();
+
+  const G4String ch = fMaterial->GetChemicalFormula();
+
+  if(ch != "") fMeanExcitationEnergy = FindMeanExcitationEnergy(ch);
+
+  // Chemical formula defines mean excitation energy
+  if(fMeanExcitationEnergy > 0.0) {
+    fLogMeanExcEnergy = std::log(fMeanExcitationEnergy);
+
+    // Compute average 
+  } else {
+    for (size_t i=0; i < nElements; i++) {
+      const G4Element* elm = (*elmVector)[i];
+      fLogMeanExcEnergy += nAtomsPerVolume[i]*elm->GetZ()
+	*std::log(elm->GetIonisation()->GetMeanExcitationEnergy());
+    }
+    fLogMeanExcEnergy /= fMaterial->GetTotNbOfElectPerVolume();
+    fMeanExcitationEnergy = std::exp(fLogMeanExcEnergy);
   }
-
-  fLogMeanExcEnergy /= fMaterial->GetTotNbOfElectPerVolume();
-  fMeanExcitationEnergy = std::exp(fLogMeanExcEnergy);
 
   fShellCorrectionVector = new G4double[3]; //[3]
@@ -96,8 +110,7 @@
     fShellCorrectionVector[j] = 0.;
 
-    for (size_t k=0; k<fMaterial->GetNumberOfElements(); k++) {
-      fShellCorrectionVector[j] += (fMaterial->GetVecNbOfAtomsPerVolume())[k] 
-              *((*(fMaterial->GetElementVector()))[k]->GetIonisation()
-                                              ->GetShellCorrectionVector()[j]);
+    for (size_t k=0; k<nElements; k++) {
+      fShellCorrectionVector[j] += nAtomsPerVolume[k]
+	*(((*elmVector)[k])->GetIonisation()->GetShellCorrectionVector())[j];
     }
     fShellCorrectionVector[j] *= 2.0/fMaterial->GetTotNbOfElectPerVolume();
@@ -200,8 +213,8 @@
   // need an 'effective Z' ?????
   G4double Zeff = 0.;
-  for (size_t i=0;i<fMaterial->GetNumberOfElements();i++) 
+  for (size_t i=0;i<fMaterial->GetNumberOfElements();i++) {
      Zeff += (fMaterial->GetFractionVector())[i]
              *((*(fMaterial->GetElementVector()))[i]->GetZ());
-
+  }
   if (Zeff > 2.) fF2fluct = 2./Zeff ;
     else         fF2fluct = 0.;
@@ -221,35 +234,4 @@
 void G4IonisParamMat::ComputeIonParameters()
 {
-  // compute parameters for ion transport
-  // The aproximation from:
-  // J.F.Ziegler, J.P. Biersack, U. Littmark
-  // The Stopping and Range of Ions in Matter,
-  // Vol.1, Pergamon Press, 1985
-  // Fast ions or hadrons
-
-  static G4double vFermi[92] = {
-    1.0309,  0.15976, 0.59782, 1.0781,  1.0486,  1.0,     1.058,   0.93942, 0.74562, 0.3424,
-    0.45259, 0.71074, 0.90519, 0.97411, 0.97184, 0.89852, 0.70827, 0.39816, 0.36552, 0.62712,
-    0.81707, 0.9943,  1.1423,  1.2381,  1.1222,  0.92705, 1.0047,  1.2,     1.0661,  0.97411,
-    0.84912, 0.95,    1.0903,  1.0429,  0.49715, 0.37755, 0.35211, 0.57801, 0.77773, 1.0207,
-    1.029,   1.2542,  1.122,   1.1241,  1.0882,  1.2709,  1.2542,  0.90094, 0.74093, 0.86054,
-    0.93155, 1.0047,  0.55379, 0.43289, 0.32636, 0.5131,  0.695,   0.72591, 0.71202, 0.67413,
-    0.71418, 0.71453, 0.5911,  0.70263, 0.68049, 0.68203, 0.68121, 0.68532, 0.68715, 0.61884,
-    0.71801, 0.83048, 1.1222,  1.2381,  1.045,   1.0733,  1.0953,  1.2381,  1.2879,  0.78654,
-    0.66401, 0.84912, 0.88433, 0.80746, 0.43357, 0.41923, 0.43638, 0.51464, 0.73087, 0.81065,
-    1.9578,  1.0257} ;
-
-  static G4double lFactor[92] = {
-    1.0,  1.0,  1.1,  1.06, 1.01, 1.03, 1.04, 0.99, 0.95, 0.9,
-    0.82, 0.81, 0.83, 0.88, 1.0,  0.95, 0.97, 0.99, 0.98, 0.97,
-    0.98, 0.97, 0.96, 0.93, 0.91, 0.9,  0.88, 0.9,  0.9,  0.9,
-    0.9,  0.85, 0.9,  0.9,  0.91, 0.92, 0.9,  0.9,  0.9,  0.9,
-    0.9,  0.88, 0.9,  0.88, 0.88, 0.9,  0.9,  0.88, 0.9,  0.9,
-    0.9,  0.9,  0.96, 1.2,  0.9,  0.88, 0.88, 0.85, 0.9,  0.9,
-    0.92, 0.95, 0.99, 1.03, 1.05, 1.07, 1.08, 1.1,  1.08, 1.08,
-    1.08, 1.08, 1.09, 1.09, 1.1,  1.11, 1.12, 1.13, 1.14, 1.15,
-    1.17, 1.2,  1.18, 1.17, 1.17, 1.16, 1.16, 1.16, 1.16, 1.16,
-    1.16, 1.16} ;
-
   // get elements in the actual material,
   const G4ElementVector* theElementVector = fMaterial->GetElementVector() ;
@@ -263,10 +245,8 @@
 
   if( 1 == NumberOfElements ) {
-    z = fMaterial->GetZ() ;
-    G4int iz = G4int(z) - 1 ;
-    if(iz < 0) iz = 0 ;
-    else if(iz > 91) iz = 91 ;
-    vF   = vFermi[iz] ;
-    lF   = lFactor[iz] ;
+    const G4Element* element = (*theElementVector)[0];
+    z = element->GetZ();
+    vF= element->GetIonisation()->GetFermiVelocity();
+    lF= element->GetIonisation()->GetLFactor();
 
   } else {
@@ -274,13 +254,9 @@
       {
         const G4Element* element = (*theElementVector)[iel] ;
-        G4double z2 = element->GetZ() ;
         const G4double weight = theAtomicNumDensityVector[iel] ;
         norm += weight ;
-        z    += z2 * weight ;
-        G4int iz = G4int(z2) - 1 ;
-        if(iz < 0) iz = 0 ;
-        else if(iz > 91) iz =91 ;
-        vF   += vFermi[iz] * weight ;
-        lF   += lFactor[iz] * weight ;
+        z    += element->GetZ() * weight ;
+        vF   += element->GetIonisation()->GetFermiVelocity() * weight ;
+        lF   += element->GetIonisation()->GetLFactor() * weight ;
       }
     z  /= norm ;
@@ -299,4 +275,5 @@
   if(value == fMeanExcitationEnergy || value <= 0.0) return;
 
+  /*
   if (G4NistManager::Instance()->GetVerbose() > 0) 
     G4cout << "G4Material: Mean excitation energy is changed for "
@@ -305,5 +282,6 @@
            << "eV; Inew= " << value/eV << " eV;"
            << G4endl;
-
+  */
+  
   fMeanExcitationEnergy = value;
   fLogMeanExcEnergy = std::log(value);
Index: trunk/source/materials/src/G4Isotope.cc
===================================================================
--- trunk/source/materials/src/G4Isotope.cc	(revision 822)
+++ trunk/source/materials/src/G4Isotope.cc	(revision 850)
@@ -25,6 +25,6 @@
 //
 //
-// $Id: G4Isotope.cc,v 1.21 2007/10/18 11:14:33 vnivanch Exp $
-// GEANT4 tag $Name:  $
+// $Id: G4Isotope.cc,v 1.22 2008/08/11 11:53:11 vnivanch Exp $
+// GEANT4 tag $Name: HEAD $
 //
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
@@ -63,7 +63,7 @@
     (" ERROR! Attempt to create an Isotope with N < Z !!!" );
     
-  if (A<=DBL_MIN)
-    fA = (G4NistManager::Instance()->GetIsotopeMass(Z,N))*g/(mole*amu_c2);  
-
+  if (A<=DBL_MIN) {
+    fA = (G4NistManager::Instance()->GetAtomicMass(Z,N))*g/(mole*amu_c2);  
+  }
   theIsotopeTable.push_back(this);
   fIndexInTable = theIsotopeTable.size() - 1;
Index: trunk/source/materials/src/G4MPVEntry.cc
===================================================================
--- trunk/source/materials/src/G4MPVEntry.cc	(revision 822)
+++ trunk/source/materials/src/G4MPVEntry.cc	(revision 850)
@@ -25,6 +25,6 @@
 //
 //
-// $Id: G4MPVEntry.cc,v 1.7 2006/06/29 19:12:49 gunter Exp $
-// GEANT4 tag $Name:  $
+// $Id: G4MPVEntry.cc,v 1.8 2008/06/05 23:39:18 gum Exp $
+// GEANT4 tag $Name: HEAD $
 //
 // 
@@ -59,5 +59,5 @@
 G4bool G4MPVEntry::operator ==(const G4MPVEntry &right) const  
 {
-	if (thePhotonMomentum == right.thePhotonMomentum) 
+	if (thePhotonEnergy == right.thePhotonEnergy) 
 		return true;
 	else
@@ -71,5 +71,5 @@
 G4bool G4MPVEntry::operator <(const G4MPVEntry &right) const  
 {
-	if (thePhotonMomentum < right.thePhotonMomentum) 
+	if (thePhotonEnergy < right.thePhotonEnergy) 
 		return true;
 	else
@@ -85,5 +85,5 @@
         if (this == &right) return *this;
 	
-        thePhotonMomentum = right.thePhotonMomentum;
+        thePhotonEnergy = right.thePhotonEnergy;
         theProperty = right.theProperty;
         return *this;
@@ -94,7 +94,7 @@
         /////////////////
 
-G4MPVEntry::G4MPVEntry(G4double aPhotonMomentum, G4double aProperty)
+G4MPVEntry::G4MPVEntry(G4double aPhotonEnergy, G4double aProperty)
 {
-        thePhotonMomentum = aPhotonMomentum;
+        thePhotonEnergy = aPhotonEnergy;
         theProperty = aProperty;
 }
@@ -102,5 +102,5 @@
 G4MPVEntry::G4MPVEntry(const G4MPVEntry &right)
 {
-        thePhotonMomentum = right.thePhotonMomentum;
+        thePhotonEnergy = right.thePhotonEnergy;
         theProperty = right.theProperty;
 }
@@ -120,5 +120,5 @@
 {
 	G4cout << "(" 
-	     << thePhotonMomentum 
+	     << thePhotonEnergy
 	     << ", " 
 	     << theProperty
Index: trunk/source/materials/src/G4Material.cc
===================================================================
--- trunk/source/materials/src/G4Material.cc	(revision 822)
+++ trunk/source/materials/src/G4Material.cc	(revision 850)
@@ -25,6 +25,6 @@
 //
 //
-// $Id: G4Material.cc,v 1.38.2.1 2008/04/28 07:24:36 gcosmo Exp $
-// GEANT4 tag $Name: geant4-09-01-patch-02 $
+// $Id: G4Material.cc,v 1.42 2008/08/13 16:06:42 vnivanch Exp $
+// GEANT4 tag $Name: HEAD $
 //
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
@@ -66,4 +66,5 @@
 // 27-07-07, improve destructor (V.Ivanchenko) 
 // 18-10-07, move definition of material index to InitialisePointers (V.Ivanchenko) 
+// 13-08-08, do not use fixed size arrays (V.Ivanchenko) 
 // 
 
@@ -107,4 +108,5 @@
   // element corresponding to this material
   maxNbComponents        = fNumberOfComponents = fNumberOfElements = 1;
+  fArrayLength           = maxNbComponents;
   fImplicitElement       = true;
   theElementVector       = new G4ElementVector();
@@ -152,4 +154,5 @@
     
   maxNbComponents     = nComponents;
+  fArrayLength        = maxNbComponents;
   fNumberOfComponents = fNumberOfElements = 0;
   fImplicitElement    = false;
@@ -185,6 +188,6 @@
   // initialization
   if ( fNumberOfElements == 0 ) {
-     fAtomsVector        = new G4int   [maxNbComponents];
-     fMassFractionVector = new G4double[maxNbComponents];
+     fAtomsVector        = new G4int   [fArrayLength];
+     fMassFractionVector = new G4double[fArrayLength];
   }
 
@@ -195,9 +198,10 @@
      fNumberOfComponents = ++fNumberOfElements;
      element->increaseCountUse();
-  }
-  else
-     G4Exception
+  } else {
+    G4cerr << "G4Material::AddElement ERROR for " << fName << " nElement= " 
+	   <<  fNumberOfElements << G4endl;
+    G4Exception
     ("ERROR!!! - Attempt to add more than the declared number of elements.");
-
+  }
   // filled.
   if ( G4int(fNumberOfElements) == maxNbComponents ) {     
@@ -226,6 +230,6 @@
   // initialization
   if (fNumberOfComponents == 0) {
-    fMassFractionVector = new G4double[50];
-    fAtomsVector        = new G4int   [50];
+    fMassFractionVector = new G4double[fArrayLength];
+    fAtomsVector        = new G4int   [fArrayLength];
   }
 
@@ -242,9 +246,11 @@
       }
       fNumberOfComponents++;  
-  }
-  else
-     G4Exception
+  } else {
+    G4cerr << "G4Material::AddElement ERROR for " << fName << " nElement= " 
+	   <<  fNumberOfElements << G4endl;
+    G4Exception
     ("ERROR!!! - Attempt to add more than the declared number of components.");
-    
+  }    
+
   // filled.
   if (G4int(fNumberOfComponents) == maxNbComponents) {
@@ -281,11 +287,29 @@
   // initialization
   if (fNumberOfComponents == 0) {
-    fMassFractionVector = new G4double[50];
-    fAtomsVector        = new G4int   [50];
+    fMassFractionVector = new G4double[fArrayLength];
+    fAtomsVector        = new G4int   [fArrayLength];
+  }
+
+  size_t nelm = material->GetNumberOfElements();
+
+  // arrays should be extended
+  if(nelm > 1) {
+    G4int nold    = fArrayLength;
+    fArrayLength += nelm - 1;
+    G4double* v1 = new G4double[fArrayLength];
+    G4int* i1    = new G4int[fArrayLength];
+    for(G4int i=0; i<nold; i++) {
+      v1[i] = fMassFractionVector[i];
+      i1[i] = fAtomsVector[i];
+    }
+    delete [] fAtomsVector;
+    delete [] fMassFractionVector;
+    fMassFractionVector = v1;
+    fAtomsVector = i1;
   }
 
   // filling ...
   if (G4int(fNumberOfComponents) < maxNbComponents) {
-     for (size_t elm=0; elm < material->GetNumberOfElements(); elm++)
+     for (size_t elm=0; elm<nelm; elm++)
        {
         G4Element* element = (*(material->GetElementVector()))[elm];
@@ -303,9 +327,11 @@
        } 
       fNumberOfComponents++;  
-  }
-  else
-     G4Exception
-    ("ERROR!!! - Attempt to add more than the declared number of components.");
-    
+  } else {
+    G4cerr << "G4Material::AddElement ERROR for " << fName << " nElement= " 
+	   <<  fNumberOfElements << G4endl;
+    G4Exception
+      ("ERROR!!! - Attempt to add more than the declared number of components.");
+  }    
+
   // filled.
   if (G4int(fNumberOfComponents) == maxNbComponents) {
@@ -541,12 +567,15 @@
   
   flux
-    << " Material: "      << std::setw(8) <<  material->fName
+    << " Material: "         << std::setw(8) <<  material->fName
     << " " << material->fChemicalFormula << " "
-    << "  density: "     << std::setw(6) << std::setprecision(3)  
+    << "  density: "         << std::setw(6) << std::setprecision(3)  
     << G4BestUnit(material->fDensity,"Volumic Mass") 
-    << "  RadL: "        << std::setw(7)  << std::setprecision(3)  
+    << "  RadL: "            << std::setw(7)  << std::setprecision(3)  
     << G4BestUnit(material->fRadlen,"Length")
-    << "  Imean: "       << std::setw(7)  << std::setprecision(3)  
+    << "  Nucl.Int.Length: " << std::setw(7)  << std::setprecision(3)  
+    << G4BestUnit(material->fNuclInterLen,"Length")    
+    << "  Imean: "           << std::setw(7)  << std::setprecision(3)  
     << G4BestUnit(material->GetIonisation()->GetMeanExcitationEnergy(),"Energy");
+    
   if(material->fState == kStateGas)
     flux
Index: trunk/source/materials/src/G4MaterialPropertiesTable.cc
===================================================================
--- trunk/source/materials/src/G4MaterialPropertiesTable.cc	(revision 822)
+++ trunk/source/materials/src/G4MaterialPropertiesTable.cc	(revision 850)
@@ -25,6 +25,6 @@
 //
 //
-// $Id: G4MaterialPropertiesTable.cc,v 1.19 2007/04/25 15:46:55 gum Exp $
-// GEANT4 tag $Name:  $
+// $Id: G4MaterialPropertiesTable.cc,v 1.20 2008/06/05 23:38:34 gum Exp $
+// GEANT4 tag $Name: HEAD $
 //
 // 
@@ -86,5 +86,5 @@
 
 void G4MaterialPropertiesTable::AddProperty(const char     *key,
-					    G4double *PhotonMomenta,
+					    G4double *PhotonEnergies,
 					    G4double *PropertyValues,
 					    G4int     NumEntries)
@@ -94,5 +94,5 @@
 
 	G4MaterialPropertyVector *mpv = 
-			new G4MaterialPropertyVector(PhotonMomenta, 
+			new G4MaterialPropertyVector(PhotonEnergies, 
 					  	     PropertyValues, 
 						     NumEntries);
@@ -165,5 +165,5 @@
 
 void G4MaterialPropertiesTable::AddEntry(const char     *key,
-					 G4double  aPhotonMomentum,
+					 G4double  aPhotonEnergy,
 					 G4double  aPropertyValue)
 
@@ -174,5 +174,5 @@
 	G4MaterialPropertyVector *targetVector=MPT [G4String(key)];
 	if (targetVector != 0) {
-		targetVector->AddElement(aPhotonMomentum, aPropertyValue);
+		targetVector->AddElement(aPhotonEnergy, aPropertyValue);
 	}
 	else {
@@ -183,5 +183,5 @@
 
 void G4MaterialPropertiesTable::RemoveEntry(const char *key,  
-					    G4double  aPhotonMomentum)
+					    G4double  aPhotonEnergy)
 {
 //      Allows to remove an entry pair directly from the Material Property Vector
@@ -190,5 +190,5 @@
         G4MaterialPropertyVector *targetVector=MPT [G4String(key)];
 	if (targetVector) {
-		targetVector->RemoveElement(aPhotonMomentum);
+		targetVector->RemoveElement(aPhotonEnergy);
  	}
         else {
@@ -238,5 +238,5 @@
   // fill GROUPVEL vector using RINDEX values
   // rindex built-in "iterator" was advanced to first entry above
-  G4double E0 = rindex->GetPhotonMomentum();
+  G4double E0 = rindex->GetPhotonEnergy();
   G4double n0 = rindex->GetProperty();
 
@@ -247,5 +247,5 @@
     // good, we have at least two entries in RINDEX
     // get next energy/value pair
-    G4double E1 = rindex->GetPhotonMomentum();
+    G4double E1 = rindex->GetPhotonEnergy();
     G4double n1 = rindex->GetProperty();
 
@@ -268,5 +268,5 @@
       E0 = E1;
       n0 = n1;
-      E1 = rindex->GetPhotonMomentum();
+      E1 = rindex->GetPhotonEnergy();
       n1 = rindex->GetProperty();
 
Index: trunk/source/materials/src/G4MaterialPropertyVector.cc
===================================================================
--- trunk/source/materials/src/G4MaterialPropertyVector.cc	(revision 822)
+++ trunk/source/materials/src/G4MaterialPropertyVector.cc	(revision 850)
@@ -25,6 +25,6 @@
 //
 //
-// $Id: G4MaterialPropertyVector.cc,v 1.15 2006/06/29 19:12:56 gunter Exp $
-// GEANT4 tag $Name:  $
+// $Id: G4MaterialPropertyVector.cc,v 1.16 2008/06/05 23:38:51 gum Exp $
+// GEANT4 tag $Name: HEAD $
 //
 // 
@@ -98,5 +98,5 @@
         /////////////////
 
-G4MaterialPropertyVector::G4MaterialPropertyVector(G4double *PhotonMomenta,
+G4MaterialPropertyVector::G4MaterialPropertyVector(G4double *PhotonEnergies,
  						   G4double *PropertyValues,
 						   G4int     NumElements)
@@ -106,8 +106,8 @@
 
 	// create a vector filling it with the values
-	// from PhotonMomenta[] and PropertyValues[] 
+	// from PhotonEnergies[] and PropertyValues[] 
 
 	for(G4int i = 0; i < NumElements; i++) {
-		AddElement(PhotonMomenta[i], PropertyValues[i]);
+		AddElement(PhotonEnergies[i], PropertyValues[i]);
         }
 }
@@ -141,10 +141,10 @@
         // Methods
         ////////////
-void G4MaterialPropertyVector::AddElement(G4double aPhotonMomentum,
+void G4MaterialPropertyVector::AddElement(G4double aPhotonEnergy,
 					  G4double aPropertyValue) 
 {
 	G4MPVEntry *newElement;
 	
-	newElement = new G4MPVEntry(aPhotonMomentum, aPropertyValue);
+	newElement = new G4MPVEntry(aPhotonEnergy, aPropertyValue);
 	MPV.push_back(newElement);
 	std::sort(MPV.begin(), MPV.end(), MPVEntry_less());
@@ -152,10 +152,10 @@
 }
 
-void G4MaterialPropertyVector::RemoveElement(G4double aPhotonMomentum)
+void G4MaterialPropertyVector::RemoveElement(G4double aPhotonEnergy)
 {
 	G4MPVEntry *newElement;
 	G4MPVEntry *success=0;
 
-	newElement = new G4MPVEntry(aPhotonMomentum, DBL_MAX);
+	newElement = new G4MPVEntry(aPhotonEnergy, DBL_MAX);
 
 	std::vector<G4MPVEntry*>::iterator i;
@@ -179,5 +179,5 @@
 
 G4double 
-G4MaterialPropertyVector::GetProperty(G4double aPhotonMomentum) const
+G4MaterialPropertyVector::GetProperty(G4double aPhotonEnergy) const
 {
 	G4MPVEntry *target, *temp; 
@@ -189,7 +189,7 @@
 	/////////////////////////
 
-	G4double PMmin   = MPV.front()->GetPhotonMomentum(); 
+	G4double PMmin   = MPV.front()->GetPhotonEnergy(); 
 	G4double minProp = MPV.front()->GetProperty(); 
-	G4double PMmax   = MPV.back() ->GetPhotonMomentum();
+	G4double PMmax   = MPV.back() ->GetPhotonEnergy();
 	G4double maxProp = MPV.back() ->GetProperty();
 
@@ -198,5 +198,5 @@
 	///////////////////////////////////////////
 
-	if (aPhotonMomentum < PMmin) 
+	if (aPhotonEnergy < PMmin) 
 	{
 		G4cout << "\nWarning: G4MaterialPropertyVector::GetProperty";
@@ -206,5 +206,5 @@
 	} 
 
-        if (aPhotonMomentum > PMmax)
+        if (aPhotonEnergy > PMmax)
 	{
 		G4cout << "\nWarning: G4MaterialPropertyVector::GetProperty";
@@ -214,5 +214,5 @@
 	} 
 	
-	target = new G4MPVEntry(aPhotonMomentum, 0.0);
+	target = new G4MPVEntry(aPhotonEnergy, 0.0);
 
 	temp = 0;
@@ -236,9 +236,9 @@
 		//////////////////////////////
 
-		GetAdjacentBins(aPhotonMomentum, &left, &right);
-
-                pmleft = MPV[left]->GetPhotonMomentum();
-                pmright = MPV[right]->GetPhotonMomentum();
-                ratio1 = (aPhotonMomentum-pmleft)/(pmright-pmleft);
+		GetAdjacentBins(aPhotonEnergy, &left, &right);
+
+                pmleft = MPV[left]->GetPhotonEnergy();
+                pmright = MPV[right]->GetPhotonEnergy();
+                ratio1 = (aPhotonEnergy-pmleft)/(pmright-pmleft);
                 ratio2 = 1 - ratio1;
 		InterpolatedValue = MPV[left]->GetProperty()*ratio2 + 
@@ -264,29 +264,29 @@
 }
 
-G4double G4MaterialPropertyVector::GetPhotonMomentum() const
+G4double G4MaterialPropertyVector::GetPhotonEnergy() const
 {
 // 	For use with G4MaterialPropertyVector iterator
 
 	if(CurrentEntry == -1 || CurrentEntry >= NumEntries) {
-		G4Exception("G4MaterialPropertyVector::GetPhotonMomentum ==>"
-		"Iterator attempted to Retrieve Photon Momentum out of range");
+		G4Exception("G4MaterialPropertyVector::GetPhotonEnergy ==>"
+		"Iterator attempted to Retrieve Photon Energy out of range");
                 return DBL_MAX;
 	}
 	else { 
-		return MPV[CurrentEntry]->GetPhotonMomentum();
+		return MPV[CurrentEntry]->GetPhotonEnergy();
 	}
 }
 
 G4double 
-G4MaterialPropertyVector::GetPhotonMomentum(G4double aProperty) const
+G4MaterialPropertyVector::GetPhotonEnergy(G4double aProperty) const
 {
 // 				***NB*** 
-// Assumes that the property is an increasing function of photon momentum (e.g.
+// Assumes that the property is an increasing function of photon energy (e.g.
 // refraction index)
 // 				***NB*** 
 //
-// Returns the photon momentum corresponding to the property value passed in.
-// If several photon momentum values correspond to the value passed in, the
-// function returns the first photon momentum in the vector that corresponds
+// Returns the photon energy corresponding to the property value passed in.
+// If several photon energy values correspond to the value passed in, the
+// function returns the first photon energy in the vector that corresponds
 // to that value. 
 
@@ -299,7 +299,7 @@
 
 	G4double PropMin = MPV.front()->GetProperty();
-	G4double PMmin   = MPV.front()->GetPhotonMomentum();
+	G4double PMmin   = MPV.front()->GetPhotonEnergy();
 	G4double PropMax = MPV.back() ->GetProperty();
-	G4double PMmax   = MPV.back() ->GetPhotonMomentum();
+	G4double PMmax   = MPV.back() ->GetPhotonEnergy();
 
 	///////////////////////////////////////////
@@ -309,6 +309,6 @@
 	if (aProperty < PropMin) 
 	{
-	  G4cout << "\nWarning: G4MaterialPropertyVector::GetPhotonMomentum";
-	  G4cout << "\n==> attempt to Retrieve Photon Momentum out of range" 
+	  G4cout << "\nWarning: G4MaterialPropertyVector::GetPhotonEnergy";
+	  G4cout << "\n==> attempt to Retrieve Photon Energy out of range" 
 	       << G4endl; 
 	  return PMmin;
@@ -316,6 +316,6 @@
 
 	if (aProperty > PropMax) {
-          G4cout << "\nWarning: G4MaterialPropertyVector::GetPhotonMomentum";
-	  G4cout << "\n==> attempt to Retrieve Photon Momentum out of range" 
+          G4cout << "\nWarning: G4MaterialPropertyVector::GetPhotonEnergy";
+	  G4cout << "\n==> attempt to Retrieve Photon Energy out of range" 
 	       << G4endl;
 	  return PMmax;
@@ -335,5 +335,5 @@
 		if (MPV[mid]->GetProperty() == aProperty) { 
 
-			// Get first photon momentum value in vector that 
+			// Get first photon energy value in vector that 
 			// corresponds to property value  
 
@@ -343,6 +343,6 @@
 			}
 
-			InterpolatedValue = MPV[mid]->GetPhotonMomentum();
-			goto end_GetPhotonMomentum;	
+			InterpolatedValue = MPV[mid]->GetPhotonEnergy();
+			goto end_GetPhotonEnergy;	
 		}
 		if (MPV[mid]->GetProperty() < aProperty)
@@ -357,8 +357,8 @@
 	ratio1 = (aProperty - pleft) / (pright - pleft);
 	ratio2 = 1 - ratio1;
-	InterpolatedValue = MPV[left]->GetPhotonMomentum()*ratio2 +
-			    MPV[right]->GetPhotonMomentum()*ratio1;
-
-end_GetPhotonMomentum:
+	InterpolatedValue = MPV[left]->GetPhotonEnergy()*ratio2 +
+			    MPV[right]->GetPhotonEnergy()*ratio1;
+
+end_GetPhotonEnergy:
 
 	return InterpolatedValue;
@@ -392,5 +392,5 @@
 
 void 
-G4MaterialPropertyVector::GetAdjacentBins(G4double aPhotonMomentum,
+G4MaterialPropertyVector::GetAdjacentBins(G4double aPhotonEnergy,
                                           G4int *left,
 					  G4int *right) const
@@ -401,10 +401,10 @@
         *right = (MPV.size() - 1); // was .entries()
 
-        // find values in bins on either side of aPhotonMomentum
+        // find values in bins on either side of aPhotonEnergy
 
 	do {
         	mid = (*left + *right)/2;
 
-                if (MPV[mid]->GetPhotonMomentum() < aPhotonMomentum) 
+                if (MPV[mid]->GetPhotonEnergy() < aPhotonEnergy) 
 		{
                 	*left = mid;
Index: trunk/source/materials/src/G4NistElementBuilder.cc
===================================================================
--- trunk/source/materials/src/G4NistElementBuilder.cc	(revision 822)
+++ trunk/source/materials/src/G4NistElementBuilder.cc	(revision 850)
@@ -24,6 +24,6 @@
 // ********************************************************************
 //
-// $Id: G4NistElementBuilder.cc,v 1.16 2007/07/28 15:58:03 vnivanch Exp $
-// GEANT4 tag $Name:  $
+// $Id: G4NistElementBuilder.cc,v 1.22 2008/08/11 11:53:11 vnivanch Exp $
+// GEANT4 tag $Name: HEAD $
 //
 // -------------------------------------------------------------------
@@ -58,4 +58,6 @@
 #include "G4NistElementBuilder.hh"
 #include "G4Element.hh"
+#include <sstream>
+
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
@@ -77,6 +79,15 @@
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 
+G4int G4NistElementBuilder::GetZ(const G4String& name)
+{
+  G4int Z = maxNumElements;
+  do {Z--;} while( Z>0 && elmSymbol[Z] != name);
+  return Z;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
 G4Element* G4NistElementBuilder::FindOrBuildElement(const G4String& symb,
-                                                          G4bool buildIsotopes)
+						    G4bool buildIsotopes)
 {
   if(first) {
@@ -140,6 +151,13 @@
     for (G4int i=0; i<nc; i++) {
        if (relAbundance[idx + i] > 0.0) {
-         ist = new G4Isotope(elmSymbol[Z],Z, n0 + i,
-	                     massIsotopes[idx + i]*gram/mole);
+	 std::ostringstream os; 
+	 os << elmSymbol[Z] << n0 + i;
+         ist = new G4Isotope(os.str(), Z, n0 + i, 
+			     GetAtomicMass(Z, n0 + i)*g/(mole*amu_c2));
+	 /*
+	 G4cout << " Z= " << Z << " N= " << n0 + i
+		<< " miso(amu)= " <<  GetIsotopeMass(Z, n0 + i)/amu_c2
+		<< " matom(amu)= " << GetAtomicMass(Z, n0 + i)/amu_c2 << G4endl;
+	 */
 	 iso.push_back(ist);
        }
@@ -193,5 +211,5 @@
       G4cout << G4endl;
       G4cout << "          mass(amu): ";
-      for(j=0; j<nc; j++) {G4cout << massIsotopes[idx + j] << " ";}
+      for(j=0; j<nc; j++) {G4cout << GetAtomicMass(i, n0 + j) << " ";}
       G4cout << G4endl;
       G4cout << "     abanbance: ";
@@ -232,12 +250,13 @@
   G4double www;
   size_t nm = nc;
-  //  G4double delm = G4double(Z)*electron_mass_c2/amu_c2;
 
   for(size_t i=0; i<nm; i++) {
     www = 0.01*(&W)[i];
-    massIsotopes[index] = (&A)[i]; // - delm;
-    sigMass[index]      = (&sA)[i];
+    // mass of the isotope in G4 units
+    massIsotopes[index] = (&A)[i]*amu_c2 - Z*electron_mass_c2 + bindingEnergy[Z]; 
+    sigMass[index]      = (&sA)[i]*amu_c2;
     relAbundance[index] = www;
 
+    // computation of mean atomic mass of the element in atomic units
     atomicMass[Z] += www*(&A)[i];
     ww += www;
@@ -258,4 +277,13 @@
 void G4NistElementBuilder::Initialise()
 {
+  // Parameterisation from D.Lunney,J.M.Pearson,C.Thibault, 
+  // Rev.Mod.Phys. 75 (2003) 1021 
+  bindingEnergy[0] = 0.0;
+  for(G4int i=1; i<maxNumElements; i++) {
+    G4double Z = G4double(i);
+    bindingEnergy[i] = (14.4381*std::pow(Z,2.39) + 1.55468e-6*std::pow(Z,5.35))*eV;
+  }
+
+  // NIST data
   index    = 0;
    
@@ -266,4 +294,9 @@
   double HA[6] = 
   {1.00783, 2.0141, 3.01605, 4.02783, 5.03954, 6.04494};
+
+  // Garantee consistence with G4 masses
+  HA[0] = (proton_mass_c2 + electron_mass_c2 - bindingEnergy[1])/amu_c2; 
+  HA[1] = (1.875613*GeV   + electron_mass_c2 - bindingEnergy[1])/amu_c2; 
+  HA[2] = (2.80925*GeV    + electron_mass_c2 - bindingEnergy[1])/amu_c2; 
   
   double HS[6] = 
@@ -281,4 +314,8 @@
   double HeA[8] = 
   {3.01603, 4.0026, 5.01222, 6.01889, 7.02803, 8.03392, 9.04382, 10.0524};
+
+  // Garantee consistence with G4 masses
+  HeA[0] = (2.80923*GeV  + 2.0*electron_mass_c2 - bindingEnergy[2])/amu_c2; 
+  HeA[1] = (3.727417*GeV + 2.0*electron_mass_c2 - bindingEnergy[2])/amu_c2; 
   
   double HeS[8] = 
Index: trunk/source/materials/src/G4NistManager.cc
===================================================================
--- trunk/source/materials/src/G4NistManager.cc	(revision 822)
+++ trunk/source/materials/src/G4NistManager.cc	(revision 850)
@@ -24,6 +24,6 @@
 // ********************************************************************
 //
-// $Id: G4NistManager.cc,v 1.15 2007/12/11 13:32:08 gcosmo Exp $
-// GEANT4 tag $Name: geant4-09-01-patch-02 $
+// $Id: G4NistManager.cc,v 1.19 2008/07/23 14:49:31 vnivanch Exp $
+// GEANT4 tag $Name: HEAD $
 //
 // -------------------------------------------------------------------
@@ -63,6 +63,4 @@
 
 G4NistManager* G4NistManager::instance = 0;
-G4double G4NistManager::POWERZ13[256] = {0};
-G4double G4NistManager::LOGA[256] = {0};
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
@@ -89,4 +87,6 @@
   
   messenger  = new G4NistMessenger(this);  
+
+  // compute frequently used values
   for(G4int i=1; i<256; i++) {
     G4double x = G4double(i);
@@ -94,6 +94,13 @@
     LOGA[i] = std::log(x);
   }
+  for(G4int j=1; j<101; j++) {
+    G4double A = elmBuilder->GetA(j);
+    POWERA27[j] = std::pow(A,0.27);
+    LOGAZ[j]    = std::log(A);
+  }
   POWERZ13[0] = 1.0;
+  POWERA27[0] = 1.0;
   LOGA[0]     = 0.0;
+  LOGAZ[0]    = 0.0;
 }
 
@@ -176,24 +183,2 @@
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-
-G4double G4NistManager::GetZ13(G4double Z)
-{
-  G4int iz = G4int(Z);
-  G4double x = (Z - G4double(iz))/(3.0*Z);
-  if(iz > 255) iz = 255;
-  else if(iz < 0) iz = 0;
-  return POWERZ13[iz]*(1.0 + x);
-}
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-
-G4double G4NistManager::GetLOGA(G4double A)
-{
-  G4int ia = G4int(A);
-  G4double x = (A - G4double(ia))/A;
-  if(ia > 255) ia = 255;
-  else if(ia < 0) ia = 0;
-  return LOGA[ia] + x;
-}
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
Index: trunk/source/materials/src/G4NistMaterialBuilder.cc
===================================================================
--- trunk/source/materials/src/G4NistMaterialBuilder.cc	(revision 822)
+++ trunk/source/materials/src/G4NistMaterialBuilder.cc	(revision 850)
@@ -24,6 +24,6 @@
 // ********************************************************************
 //
-// $Id: G4NistMaterialBuilder.cc,v 1.17.2.3 2008/04/28 08:39:03 gcosmo Exp $
-// GEANT4 tag $Name: geant4-09-01-patch-02 $
+// $Id: G4NistMaterialBuilder.cc,v 1.19 2008/04/28 08:51:29 vnivanch Exp $
+// GEANT4 tag $Name: HEAD $
 //
 //
@@ -212,5 +212,9 @@
     return 0;
   } 
+
+  // add parameters of material into internal vectors
+  // density in g/cm3, mean ionisation potential is not defined
   AddMaterial(name,dens*cm3/g,0,0.,nm,state,temp,pressure);
+
   for (G4int i=0; i<nm; i++) {
     G4int Z = G4int((elmBuilder->FindOrBuildElement(elm[i]))->GetZ());
@@ -240,5 +244,9 @@
     return 0;
   } 
+
+  // add parameters of material into internal vectors
+  // density in g/cm3, mean ionisation potential is not defined
   AddMaterial(name,dens*cm3/g,0,0.,nm,state,temp,pressure);
+
   for (G4int i=0; i<nm; i++) {
     G4int Z = G4int((elmBuilder->FindOrBuildElement(elm[i]))->GetZ());
@@ -418,4 +426,7 @@
 					G4double temp, G4double pres)
 {
+  // add parameters of material into internal vectors
+  // density in g/cm3, mean ionisation potential in eV
+
   if (nCurrent != 0) {
     G4cout << "WARNING: G4NistMaterialBuilder::AddMaterial problem: previous "
Index: trunk/source/materials/src/G4NistMessenger.cc
===================================================================
--- trunk/source/materials/src/G4NistMessenger.cc	(revision 822)
+++ trunk/source/materials/src/G4NistMessenger.cc	(revision 850)
@@ -26,5 +26,5 @@
 //
 // $Id: G4NistMessenger.cc,v 1.4 2007/05/02 10:48:52 vnivanch Exp $
-// GEANT4 tag $Name:  $
+// GEANT4 tag $Name: HEAD $
 //
 //
Index: trunk/source/materials/src/G4OpticalSurface.cc
===================================================================
--- trunk/source/materials/src/G4OpticalSurface.cc	(revision 822)
+++ trunk/source/materials/src/G4OpticalSurface.cc	(revision 850)
@@ -25,6 +25,6 @@
 //
 //
-// $Id: G4OpticalSurface.cc,v 1.10 2006/06/29 19:13:08 gunter Exp $
-// GEANT4 tag $Name:  $
+// $Id: G4OpticalSurface.cc,v 1.11 2008/07/21 20:55:51 gum Exp $
+// GEANT4 tag $Name: HEAD $
 //
 // 
@@ -104,5 +104,8 @@
 }
 
-G4OpticalSurface::~G4OpticalSurface(){}
+G4OpticalSurface::~G4OpticalSurface()
+{
+  Overwrite();
+}
 
 G4int G4OpticalSurface::operator==(const G4OpticalSurface &right) const
Index: trunk/source/materials/src/G4SandiaTable.cc
===================================================================
--- trunk/source/materials/src/G4SandiaTable.cc	(revision 822)
+++ trunk/source/materials/src/G4SandiaTable.cc	(revision 850)
@@ -26,5 +26,5 @@
 //
 // $Id: G4SandiaTable.cc,v 1.34 2007/10/02 10:13:33 vnivanch Exp $
-// GEANT4 tag $Name:  $
+// GEANT4 tag $Name: HEAD $
 //
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
Index: trunk/source/materials/src/G4SurfaceProperty.cc
===================================================================
--- trunk/source/materials/src/G4SurfaceProperty.cc	(revision 822)
+++ trunk/source/materials/src/G4SurfaceProperty.cc	(revision 850)
@@ -25,6 +25,6 @@
 //
 //
-// $Id: G4SurfaceProperty.cc,v 1.1 2007/04/25 16:19:26 gum Exp $
-// GEANT4 tag $Name:  $
+// $Id: G4SurfaceProperty.cc,v 1.4 2008/07/21 20:55:34 gum Exp $
+// GEANT4 tag $Name: HEAD $
 //
 // 
@@ -52,4 +52,19 @@
 
 //
+// Constructor and destructor
+//
+G4SurfaceProperty::G4SurfaceProperty( const G4String& name,
+			                    G4SurfaceType type )
+  : theName(name), theType(type)
+{
+  theSurfacePropertyTable.push_back(this);
+}
+
+G4SurfaceProperty::~G4SurfaceProperty()
+{
+  Overwrite();
+}
+
+//
 // Methods
 //
@@ -67,5 +82,5 @@
 // Dump info for known surface properties
 //
-void G4SurfaceProperty::DumpInfo()
+void G4SurfaceProperty::DumpTableInfo()
 {
   G4cout << "***** Surface Property Table : Nb of Surface Properties = "
@@ -85,5 +100,5 @@
 void G4SurfaceProperty::CleanSurfacePropertyTable()
 {
-  DumpInfo();
+  DumpTableInfo();
   G4SurfacePropertyTable::iterator pos;
   for(pos=theSurfacePropertyTable.begin();
@@ -93,4 +108,4 @@
   }
   theSurfacePropertyTable.clear();
-  DumpInfo();
+  DumpTableInfo();
 }
