Ignore:
Timestamp:
Sep 10, 2008, 5:40:37 PM (16 years ago)
Author:
garnier
Message:

geant4.8.2 beta

Location:
trunk/source/materials/src
Files:
16 edited

Legend:

Unmodified
Added
Removed

  • trunk/source/materials/src/G4AtomicShells.cc

    r822 r850  
    2626//
    2727// $Id: G4AtomicShells.cc,v 1.7 2006/10/17 15:15:46 vnivanch Exp $
    28 // GEANT4 tag $Name: $
     28// GEANT4 tag $Name: HEAD $
    2929//
    3030//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....

  • trunk/source/materials/src/G4Element.cc

    r822 r850  
    2525//
    2626//
    27 // $Id: G4Element.cc,v 1.26 2007/10/18 11:14:33 vnivanch Exp $
    28 // GEANT4 tag $Name: $
     27// $Id: G4Element.cc,v 1.31 2008/08/11 11:53:11 vnivanch Exp $
     28// GEANT4 tag $Name: HEAD $
    2929//
    3030//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    … …  
    6969  : fName(name), fSymbol(symbol)                     
    7070{
    71   if (zeff<1.) G4Exception (" ERROR from G4Element::G4Element !"
    72                             " It is not allowed to create an Element with Z < 1" );
    73 
    74   if (aeff/(g/mole)<zeff) G4Exception (" ERROR from G4Element::G4Element !"
    75                                        " Attempt to create an Element with N < Z !!!" );
    76 
    77   if ((zeff-G4int(zeff)) > perMillion)
     71  G4int iz = (G4int)zeff;
     72  if (zeff<1.) {
     73    G4cout << "G4Element ERROR:  " << name << " Z= " << zeff
     74           << " A= " << aeff/(g/mole) << G4endl;
     75    G4Exception (" ERROR from G4Element::G4Element !"
     76                 " It is not allowed to create an Element with Z < 1" );
     77  }
     78  if (std::abs(zeff - iz) > perMillion) {
     79    G4cout << "G4Element Warning:  " << name << " Z= " << zeff
     80           << " A= " << aeff/(g/mole) << G4endl;
    7881    G4cerr << name << " : WARNING from G4Element::G4Element !" 
    7982      " Trying to define an element as a mixture directly via effective Z."
    8083           << G4endl;
     84  }
    8185
    8286  InitializePointers();
    … …  
    8589  fNeff   = aeff/(g/mole);
    8690  fAeff   = aeff;
     91
     92  if(fNeff < 1.0) fNeff = 1.0;
     93
     94  if (fNeff < zeff) {
     95    G4cout << "G4Element ERROR:  " << name << " Z= " << zeff
     96           << " A= " << fNeff << G4endl;
     97    G4Exception (" ERROR from G4Element::G4Element !"
     98                 " Attempt to create an Element with N < Z !!!" );
     99  }
    87100   
    88   fNbOfAtomicShells = G4AtomicShells::GetNumberOfShells((G4int)fZeff);
     101  fNbOfAtomicShells = G4AtomicShells::GetNumberOfShells(iz);
    89102  fAtomicShells     = new G4double[fNbOfAtomicShells];
    90   for (G4int i=0;i<fNbOfAtomicShells;i++)
    91     fAtomicShells[i] = G4AtomicShells::GetBindingEnergy((G4int)fZeff,i);
    92            
     103  for (G4int i=0;i<fNbOfAtomicShells;i++) {
     104    fAtomicShells[i] = G4AtomicShells::GetBindingEnergy(iz, i);
     105  }
    93106  ComputeDerivedQuantities();
    94107}
    … …  
    117130void G4Element::AddIsotope(G4Isotope* isotope, G4double abundance)
    118131{
    119   if (theIsotopeVector == 0)
     132  if (theIsotopeVector == 0) {
     133    G4cout << "G4Element ERROR:  " << fName << G4endl;
    120134    G4Exception ("ERROR from G4Element::AddIsotope!"
    121135                 " Trying to add an Isotope before contructing the element.");
     136  }
     137  G4int iz = isotope->GetZ();
    122138
    123139  // filling ...
    124140  if ( fNumberOfIsotopes < theIsotopeVector->size() ) {
    125141    // check same Z
    126     if (fNumberOfIsotopes==0) fZeff = G4double(isotope->GetZ());
    127     else if (G4double(isotope->GetZ()) != fZeff)
     142    if (fNumberOfIsotopes==0) fZeff = G4double(iz);
     143    else if (G4double(iz) != fZeff) {
    128144      G4Exception ("ERROR from G4Element::AddIsotope!"
    129145                   " Try to add isotopes with different Z");
     146    }
    130147    //Z ok   
    131148    fRelativeAbundanceVector[fNumberOfIsotopes] = abundance;
    … …  
    133150    ++fNumberOfIsotopes;
    134151    isotope->increaseCountUse();
    135   }
    136   else G4Exception ("ERROR from G4Element::AddIsotope!" 
    137                     " Attempt to add more than the declared number of isotopes.");
     152  } else {
     153    G4cout << "G4Element ERROR:  " << fName << G4endl;
     154    G4Exception ("ERROR from G4Element::AddIsotope!" 
     155                 " Attempt to add more than the declared number of isotopes.");
     156  }
    138157
    139158  // filled.
    … …  
    151170    fAeff /=  wtSum;
    152171     
    153     fNbOfAtomicShells = G4AtomicShells::GetNumberOfShells((G4int)fZeff);
     172    fNbOfAtomicShells = G4AtomicShells::GetNumberOfShells(iz);
    154173    fAtomicShells     = new G4double[fNbOfAtomicShells];
    155     for (G4int j=0;j<fNbOfAtomicShells;j++)
    156       fAtomicShells[j] = G4AtomicShells::GetBindingEnergy((G4int)fZeff,j);
     174    for (G4int j=0;j<fNbOfAtomicShells;j++) {
     175      fAtomicShells[j] = G4AtomicShells::GetBindingEnergy(iz, j);
     176    }
    157177         
    158178    ComputeDerivedQuantities();
    … …  
    213233  // check if elements with same Z already exist
    214234  fIndexZ = 0;
    215   for (size_t J=0 ; J<fIndexInTable ; J++)
     235  for (size_t J=0 ; J<fIndexInTable ; J++) {
    216236    if (theElementTable[J]->GetZ() == fZeff) fIndexZ++;
    217 
     237  }
    218238  //nb of materials which use this element
    219239  fCountUse = 0;
    … …  
    257277
    258278  G4double Lrad, Lprad;
    259   G4int iz = (int)(fZeff+0.5) - 1 ;
     279  G4int iz = (G4int)(fZeff+0.5) - 1 ;
    260280  if (iz <= 3) { Lrad = Lrad_light[iz] ;  Lprad = Lprad_light[iz] ; }
    261281    else { Lrad = std::log(184.15) - logZ3 ; Lprad = std::log(1194.) - 2*logZ3;}

  • trunk/source/materials/src/G4IonisParamElm.cc

    r822 r850  
    2525//
    2626//
    27 // $Id: G4IonisParamElm.cc,v 1.14 2006/06/29 19:12:43 gunter Exp $
    28 // GEANT4 tag $Name: $
     27// $Id: G4IonisParamElm.cc,v 1.15 2008/06/03 14:30:10 vnivanch Exp $
     28// GEANT4 tag $Name: HEAD $
    2929//
    3030//
    … …  
    4747G4IonisParamElm::G4IonisParamElm(G4double Z)
    4848{
    49     if (Z < 1.) G4Exception
    50       (" ERROR! It is not allowed to create an Element with Z < 1" );
     49  if (Z < 1.)
     50    G4Exception (" ERROR! It is not allowed to create an Element with Z < 1" );
    5151   
    52     // some basic functions of the atomic number
    53     fZ = Z;
    54     fZ3  = std::pow(fZ, 1./3.);
    55     fZZ3 = std::pow(fZ*(fZ+1.), 1./3.);
    56     flogZ3 = std::log(fZ)/3.;
     52  // some basic functions of the atomic number
     53  fZ = Z;
     54  fZ3  = std::pow(fZ, 1./3.);
     55  fZZ3 = std::pow(fZ*(fZ+1.), 1./3.);
     56  flogZ3 = std::log(fZ)/3.;
    5757   
    58      // Parameters for energy loss by ionisation   
    59 
    60      // Mean excitation energy
    61      // from "Stopping Powers for Electrons and Positrons"
    62      // ICRU Report N#37, 1984  (energy in eV)
    63     static double exc[100] = {
     58  // Parameters for energy loss by ionisation   
     59
     60  // Mean excitation energy
     61  // from "Stopping Powers for Electrons and Positrons"
     62  // ICRU Report N#37, 1984  (energy in eV)
     63  static double exc[100] = {
    6464    21.8, 20.9, 13.3, 15.9, 15.2, 13.0, 11.7, 11.9, 11.5, 13.7,
    6565    13.6, 13.0, 12.8, 12.4, 11.5, 11.3, 10.2, 10.4, 10.0,  9.6,
    … …  
    7373     9.6,  9.7,  9.7,  9.8,  9.8,  9.8,  9.8,  9.9,  9.9,  9.9 };
    7474
    75     G4int iz = (G4int)Z - 1 ;
    76     if(0  > iz) iz = 0;
    77     else if(99 < iz) iz = 99 ;
    78     fMeanExcitationEnergy = fZ * exc[iz] * eV ;
    79 
    80 
    81     fTau0 = 0.1*fZ3*MeV/proton_mass_c2;
    82     fTaul = 2.*MeV/proton_mass_c2;
    83 
    84     // compute the Bethe-Bloch formula for energy = fTaul*particle mass
    85     G4double rate = fMeanExcitationEnergy/electron_mass_c2 ;
    86     G4double w = fTaul*(fTaul+2.) ;
    87     fBetheBlochLow = (fTaul+1.)*(fTaul+1.)*std::log(2.*w/rate)/w - 1. ;
    88     fBetheBlochLow = 2.*fZ*twopi_mc2_rcl2*fBetheBlochLow ;
     75  G4int iz = (G4int)Z - 1 ;
     76  if(0  > iz) iz = 0;
     77  else if(99 < iz) iz = 99 ;
     78  fMeanExcitationEnergy = fZ * exc[iz] * eV ;
     79
     80  // compute parameters for ion transport
     81  // The aproximation from:
     82  // J.F.Ziegler, J.P. Biersack, U. Littmark
     83  // The Stopping and Range of Ions in Matter,
     84  // Vol.1, Pergamon Press, 1985
     85  // Fast ions or hadrons
     86
     87  if(91 < iz) iz = 91;
     88
     89  static G4double vFermi[92] = {
     90    1.0309,  0.15976, 0.59782, 1.0781,  1.0486,  1.0,     1.058,   0.93942, 0.74562, 0.3424,
     91    0.45259, 0.71074, 0.90519, 0.97411, 0.97184, 0.89852, 0.70827, 0.39816, 0.36552, 0.62712,
     92    0.81707, 0.9943,  1.1423,  1.2381,  1.1222,  0.92705, 1.0047,  1.2,     1.0661,  0.97411,
     93    0.84912, 0.95,    1.0903,  1.0429,  0.49715, 0.37755, 0.35211, 0.57801, 0.77773, 1.0207,
     94    1.029,   1.2542,  1.122,   1.1241,  1.0882,  1.2709,  1.2542,  0.90094, 0.74093, 0.86054,
     95    0.93155, 1.0047,  0.55379, 0.43289, 0.32636, 0.5131,  0.695,   0.72591, 0.71202, 0.67413,
     96    0.71418, 0.71453, 0.5911,  0.70263, 0.68049, 0.68203, 0.68121, 0.68532, 0.68715, 0.61884,
     97    0.71801, 0.83048, 1.1222,  1.2381,  1.045,   1.0733,  1.0953,  1.2381,  1.2879,  0.78654,
     98    0.66401, 0.84912, 0.88433, 0.80746, 0.43357, 0.41923, 0.43638, 0.51464, 0.73087, 0.81065,
     99    1.9578,  1.0257} ;
     100
     101  static G4double lFactor[92] = {
     102    1.0,  1.0,  1.1,  1.06, 1.01, 1.03, 1.04, 0.99, 0.95, 0.9,
     103    0.82, 0.81, 0.83, 0.88, 1.0,  0.95, 0.97, 0.99, 0.98, 0.97,
     104    0.98, 0.97, 0.96, 0.93, 0.91, 0.9,  0.88, 0.9,  0.9,  0.9,
     105    0.9,  0.85, 0.9,  0.9,  0.91, 0.92, 0.9,  0.9,  0.9,  0.9,
     106    0.9,  0.88, 0.9,  0.88, 0.88, 0.9,  0.9,  0.88, 0.9,  0.9,
     107    0.9,  0.9,  0.96, 1.2,  0.9,  0.88, 0.88, 0.85, 0.9,  0.9,
     108    0.92, 0.95, 0.99, 1.03, 1.05, 1.07, 1.08, 1.1,  1.08, 1.08,
     109    1.08, 1.08, 1.09, 1.09, 1.1,  1.11, 1.12, 1.13, 1.14, 1.15,
     110    1.17, 1.2,  1.18, 1.17, 1.17, 1.16, 1.16, 1.16, 1.16, 1.16,
     111    1.16, 1.16} ;
     112
     113  fVFermi  = vFermi[iz];
     114  fLFactor = lFactor[iz];
     115
     116  // obsolete parameters for ionisation
     117  fTau0 = 0.1*fZ3*MeV/proton_mass_c2;
     118  fTaul = 2.*MeV/proton_mass_c2;
     119
     120  // compute the Bethe-Bloch formula for energy = fTaul*particle mass
     121  G4double rate = fMeanExcitationEnergy/electron_mass_c2 ;
     122  G4double w = fTaul*(fTaul+2.) ;
     123  fBetheBlochLow = (fTaul+1.)*(fTaul+1.)*std::log(2.*w/rate)/w - 1. ;
     124  fBetheBlochLow = 2.*fZ*twopi_mc2_rcl2*fBetheBlochLow ;
    89125 
    90     fClow = std::sqrt(fTaul)*fBetheBlochLow;
    91     fAlow = 6.458040 * fClow/fTau0;
    92     G4double Taum = 0.035*fZ3*MeV/proton_mass_c2;
    93     fBlow =-3.229020*fClow/(fTau0*std::sqrt(Taum));
    94 
    95     // Shell correction factors
    96     fShellCorrectionVector = new G4double[3]; //[3]
    97     rate = 0.001*fMeanExcitationEnergy/eV;
    98     G4double rate2 = rate*rate;
     126  fClow = std::sqrt(fTaul)*fBetheBlochLow;
     127  fAlow = 6.458040 * fClow/fTau0;
     128  G4double Taum = 0.035*fZ3*MeV/proton_mass_c2;
     129  fBlow =-3.229020*fClow/(fTau0*std::sqrt(Taum));
     130
     131  // Shell correction parameterization
     132  fShellCorrectionVector = new G4double[3]; //[3]
     133  rate = 0.001*fMeanExcitationEnergy/eV;
     134  G4double rate2 = rate*rate;
    99135    /*   
    100136    fShellCorrectionVector[0] = ( 1.10289e5 + 5.14781e8*rate)*rate2 ;
    … …  
    102138    fShellCorrectionVector[2] = (-9.92256e1 + 2.10823e5*rate)*rate2 ;
    103139    */
    104     fShellCorrectionVector[0] = ( 0.422377   + 3.858019*rate)*rate2 ;
    105     fShellCorrectionVector[1] = ( 0.0304043  - 0.1667989*rate)*rate2 ;
    106     fShellCorrectionVector[2] = (-0.00038106 + 0.00157955*rate)*rate2 ;
     140  fShellCorrectionVector[0] = ( 0.422377   + 3.858019*rate)*rate2 ;
     141  fShellCorrectionVector[1] = ( 0.0304043  - 0.1667989*rate)*rate2 ;
     142  fShellCorrectionVector[2] = (-0.00038106 + 0.00157955*rate)*rate2 ;
    107143}
    108144

  • trunk/source/materials/src/G4IonisParamMat.cc

    r822 r850  
    2525//
    2626//
    27 // $Id: G4IonisParamMat.cc,v 1.20 2007/09/27 14:05:47 vnivanch Exp $
    28 // GEANT4 tag $Name: geant4-09-01-patch-02 $
     27// $Id: G4IonisParamMat.cc,v 1.25 2008/07/08 10:34:56 vnivanch Exp $
     28// GEANT4 tag $Name: HEAD $
    2929//
    3030//
    … …  
    3939// 10-05-05, add a missing coma in FindMeanExcitationEnergy() - Bug#746 (mma)
    4040// 27-09-07, add computation of parameters for ions (V.Ivanchenko)
     41// 04-03-08, remove reference to G4NistManager. Add fBirks constant (mma)
    4142
    4243//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
    … …  
    4445#include "G4IonisParamMat.hh"
    4546#include "G4Material.hh"
    46 #include "G4NistManager.hh"
    4747
    4848//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
    … …  
    5555  ComputeFluctModel();
    5656  ComputeIonParameters();
     57 
     58  fBirks = 0.;
     59  fMeanEnergyPerIon = 0.;
    5760}
    5861
    … …  
    7881  fLogMeanExcEnergy = 0.;
    7982
    80 
    81   for (size_t i=0; i < fMaterial->GetNumberOfElements(); i++) {
    82     fLogMeanExcEnergy +=
    83              (fMaterial->GetVecNbOfAtomsPerVolume())[i]
    84             *((*(fMaterial->GetElementVector()))[i]->GetZ())
    85             *std::log((*(fMaterial->GetElementVector()))[i]->GetIonisation()
    86              ->GetMeanExcitationEnergy());
     83  size_t nElements = fMaterial->GetNumberOfElements();
     84  const G4ElementVector* elmVector = fMaterial->GetElementVector();
     85  const G4double* nAtomsPerVolume = fMaterial->GetVecNbOfAtomsPerVolume();
     86
     87  const G4String ch = fMaterial->GetChemicalFormula();
     88
     89  if(ch != "") fMeanExcitationEnergy = FindMeanExcitationEnergy(ch);
     90
     91  // Chemical formula defines mean excitation energy
     92  if(fMeanExcitationEnergy > 0.0) {
     93    fLogMeanExcEnergy = std::log(fMeanExcitationEnergy);
     94
     95    // Compute average
     96  } else {
     97    for (size_t i=0; i < nElements; i++) {
     98      const G4Element* elm = (*elmVector)[i];
     99      fLogMeanExcEnergy += nAtomsPerVolume[i]*elm->GetZ()
     100        *std::log(elm->GetIonisation()->GetMeanExcitationEnergy());
     101    }
     102    fLogMeanExcEnergy /= fMaterial->GetTotNbOfElectPerVolume();
     103    fMeanExcitationEnergy = std::exp(fLogMeanExcEnergy);
    87104  }
    88 
    89   fLogMeanExcEnergy /= fMaterial->GetTotNbOfElectPerVolume();
    90   fMeanExcitationEnergy = std::exp(fLogMeanExcEnergy);
    91105
    92106  fShellCorrectionVector = new G4double[3]; //[3]
    … …  
    96110    fShellCorrectionVector[j] = 0.;
    97111
    98     for (size_t k=0; k<fMaterial->GetNumberOfElements(); k++) {
    99       fShellCorrectionVector[j] += (fMaterial->GetVecNbOfAtomsPerVolume())[k]
    100               *((*(fMaterial->GetElementVector()))[k]->GetIonisation()
    101                                               ->GetShellCorrectionVector()[j]);
     112    for (size_t k=0; k<nElements; k++) {
     113      fShellCorrectionVector[j] += nAtomsPerVolume[k]
     114        *(((*elmVector)[k])->GetIonisation()->GetShellCorrectionVector())[j];
    102115    }
    103116    fShellCorrectionVector[j] *= 2.0/fMaterial->GetTotNbOfElectPerVolume();
    … …  
    200213  // need an 'effective Z' ?????
    201214  G4double Zeff = 0.;
    202   for (size_t i=0;i<fMaterial->GetNumberOfElements();i++)
     215  for (size_t i=0;i<fMaterial->GetNumberOfElements();i++) {
    203216     Zeff += (fMaterial->GetFractionVector())[i]
    204217             *((*(fMaterial->GetElementVector()))[i]->GetZ());
    205 
     218  }
    206219  if (Zeff > 2.) fF2fluct = 2./Zeff ;
    207220    else         fF2fluct = 0.;
    … …  
    221234void G4IonisParamMat::ComputeIonParameters()
    222235{
    223   // compute parameters for ion transport
    224   // The aproximation from:
    225   // J.F.Ziegler, J.P. Biersack, U. Littmark
    226   // The Stopping and Range of Ions in Matter,
    227   // Vol.1, Pergamon Press, 1985
    228   // Fast ions or hadrons
    229 
    230   static G4double vFermi[92] = {
    231     1.0309,  0.15976, 0.59782, 1.0781,  1.0486,  1.0,     1.058,   0.93942, 0.74562, 0.3424,
    232     0.45259, 0.71074, 0.90519, 0.97411, 0.97184, 0.89852, 0.70827, 0.39816, 0.36552, 0.62712,
    233     0.81707, 0.9943,  1.1423,  1.2381,  1.1222,  0.92705, 1.0047,  1.2,     1.0661,  0.97411,
    234     0.84912, 0.95,    1.0903,  1.0429,  0.49715, 0.37755, 0.35211, 0.57801, 0.77773, 1.0207,
    235     1.029,   1.2542,  1.122,   1.1241,  1.0882,  1.2709,  1.2542,  0.90094, 0.74093, 0.86054,
    236     0.93155, 1.0047,  0.55379, 0.43289, 0.32636, 0.5131,  0.695,   0.72591, 0.71202, 0.67413,
    237     0.71418, 0.71453, 0.5911,  0.70263, 0.68049, 0.68203, 0.68121, 0.68532, 0.68715, 0.61884,
    238     0.71801, 0.83048, 1.1222,  1.2381,  1.045,   1.0733,  1.0953,  1.2381,  1.2879,  0.78654,
    239     0.66401, 0.84912, 0.88433, 0.80746, 0.43357, 0.41923, 0.43638, 0.51464, 0.73087, 0.81065,
    240     1.9578,  1.0257} ;
    241 
    242   static G4double lFactor[92] = {
    243     1.0,  1.0,  1.1,  1.06, 1.01, 1.03, 1.04, 0.99, 0.95, 0.9,
    244     0.82, 0.81, 0.83, 0.88, 1.0,  0.95, 0.97, 0.99, 0.98, 0.97,
    245     0.98, 0.97, 0.96, 0.93, 0.91, 0.9,  0.88, 0.9,  0.9,  0.9,
    246     0.9,  0.85, 0.9,  0.9,  0.91, 0.92, 0.9,  0.9,  0.9,  0.9,
    247     0.9,  0.88, 0.9,  0.88, 0.88, 0.9,  0.9,  0.88, 0.9,  0.9,
    248     0.9,  0.9,  0.96, 1.2,  0.9,  0.88, 0.88, 0.85, 0.9,  0.9,
    249     0.92, 0.95, 0.99, 1.03, 1.05, 1.07, 1.08, 1.1,  1.08, 1.08,
    250     1.08, 1.08, 1.09, 1.09, 1.1,  1.11, 1.12, 1.13, 1.14, 1.15,
    251     1.17, 1.2,  1.18, 1.17, 1.17, 1.16, 1.16, 1.16, 1.16, 1.16,
    252     1.16, 1.16} ;
    253 
    254236  // get elements in the actual material,
    255237  const G4ElementVector* theElementVector = fMaterial->GetElementVector() ;
    … …  
    263245
    264246  if( 1 == NumberOfElements ) {
    265     z = fMaterial->GetZ() ;
    266     G4int iz = G4int(z) - 1 ;
    267     if(iz < 0) iz = 0 ;
    268     else if(iz > 91) iz = 91 ;
    269     vF   = vFermi[iz] ;
    270     lF   = lFactor[iz] ;
     247    const G4Element* element = (*theElementVector)[0];
     248    z = element->GetZ();
     249    vF= element->GetIonisation()->GetFermiVelocity();
     250    lF= element->GetIonisation()->GetLFactor();
    271251
    272252  } else {
    … …  
    274254      {
    275255        const G4Element* element = (*theElementVector)[iel] ;
    276         G4double z2 = element->GetZ() ;
    277256        const G4double weight = theAtomicNumDensityVector[iel] ;
    278257        norm += weight ;
    279         z    += z2 * weight ;
    280         G4int iz = G4int(z2) - 1 ;
    281         if(iz < 0) iz = 0 ;
    282         else if(iz > 91) iz =91 ;
    283         vF   += vFermi[iz] * weight ;
    284         lF   += lFactor[iz] * weight ;
     258        z    += element->GetZ() * weight ;
     259        vF   += element->GetIonisation()->GetFermiVelocity() * weight ;
     260        lF   += element->GetIonisation()->GetLFactor() * weight ;
    285261      }
    286262    z  /= norm ;
    … …  
    299275  if(value == fMeanExcitationEnergy || value <= 0.0) return;
    300276
     277  /*
    301278  if (G4NistManager::Instance()->GetVerbose() > 0)
    302279    G4cout << "G4Material: Mean excitation energy is changed for "
    … …  
    305282           << "eV; Inew= " << value/eV << " eV;"
    306283           << G4endl;
    307 
     284  */
     285 
    308286  fMeanExcitationEnergy = value;
    309287  fLogMeanExcEnergy = std::log(value);

  • trunk/source/materials/src/G4Isotope.cc

    r822 r850  
    2525//
    2626//
    27 // $Id: G4Isotope.cc,v 1.21 2007/10/18 11:14:33 vnivanch Exp $
    28 // GEANT4 tag $Name: $
     27// $Id: G4Isotope.cc,v 1.22 2008/08/11 11:53:11 vnivanch Exp $
     28// GEANT4 tag $Name: HEAD $
    2929//
    3030//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    … …  
    6363    (" ERROR! Attempt to create an Isotope with N < Z !!!" );
    6464   
    65   if (A<=DBL_MIN)
    66     fA = (G4NistManager::Instance()->GetIsotopeMass(Z,N))*g/(mole*amu_c2); 
    67 
     65  if (A<=DBL_MIN) {
     66    fA = (G4NistManager::Instance()->GetAtomicMass(Z,N))*g/(mole*amu_c2); 
     67  }
    6868  theIsotopeTable.push_back(this);
    6969  fIndexInTable = theIsotopeTable.size() - 1;

  • trunk/source/materials/src/G4MPVEntry.cc

    r822 r850  
    2525//
    2626//
    27 // $Id: G4MPVEntry.cc,v 1.7 2006/06/29 19:12:49 gunter Exp $
    28 // GEANT4 tag $Name: $
     27// $Id: G4MPVEntry.cc,v 1.8 2008/06/05 23:39:18 gum Exp $
     28// GEANT4 tag $Name: HEAD $
    2929//
    3030//
    … …  
    5959G4bool G4MPVEntry::operator ==(const G4MPVEntry &right) const 
    6060{
    61         if (thePhotonMomentum == right.thePhotonMomentum)
     61        if (thePhotonEnergy == right.thePhotonEnergy)
    6262                return true;
    6363        else
    … …  
    7171G4bool G4MPVEntry::operator <(const G4MPVEntry &right) const 
    7272{
    73         if (thePhotonMomentum < right.thePhotonMomentum)
     73        if (thePhotonEnergy < right.thePhotonEnergy)
    7474                return true;
    7575        else
    … …  
    8585        if (this == &right) return *this;
    8686       
    87         thePhotonMomentum = right.thePhotonMomentum;
     87        thePhotonEnergy = right.thePhotonEnergy;
    8888        theProperty = right.theProperty;
    8989        return *this;
    … …  
    9494        /////////////////
    9595
    96 G4MPVEntry::G4MPVEntry(G4double aPhotonMomentum, G4double aProperty)
     96G4MPVEntry::G4MPVEntry(G4double aPhotonEnergy, G4double aProperty)
    9797{
    98         thePhotonMomentum = aPhotonMomentum;
     98        thePhotonEnergy = aPhotonEnergy;
    9999        theProperty = aProperty;
    100100}
    … …  
    102102G4MPVEntry::G4MPVEntry(const G4MPVEntry &right)
    103103{
    104         thePhotonMomentum = right.thePhotonMomentum;
     104        thePhotonEnergy = right.thePhotonEnergy;
    105105        theProperty = right.theProperty;
    106106}
    … …  
    120120{
    121121        G4cout << "("
    122              << thePhotonMomentum
     122             << thePhotonEnergy
    123123             << ", "
    124124             << theProperty

  • trunk/source/materials/src/G4Material.cc

    r822 r850  
    2525//
    2626//
    27 // $Id: G4Material.cc,v 1.38.2.1 2008/04/28 07:24:36 gcosmo Exp $
    28 // GEANT4 tag $Name: geant4-09-01-patch-02 $
     27// $Id: G4Material.cc,v 1.42 2008/08/13 16:06:42 vnivanch Exp $
     28// GEANT4 tag $Name: HEAD $
    2929//
    3030//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    … …  
    6666// 27-07-07, improve destructor (V.Ivanchenko)
    6767// 18-10-07, move definition of material index to InitialisePointers (V.Ivanchenko)
     68// 13-08-08, do not use fixed size arrays (V.Ivanchenko)
    6869//
    6970
    … …  
    107108  // element corresponding to this material
    108109  maxNbComponents        = fNumberOfComponents = fNumberOfElements = 1;
     110  fArrayLength           = maxNbComponents;
    109111  fImplicitElement       = true;
    110112  theElementVector       = new G4ElementVector();
    … …  
    152154   
    153155  maxNbComponents     = nComponents;
     156  fArrayLength        = maxNbComponents;
    154157  fNumberOfComponents = fNumberOfElements = 0;
    155158  fImplicitElement    = false;
    … …  
    185188  // initialization
    186189  if ( fNumberOfElements == 0 ) {
    187      fAtomsVector        = new G4int   [maxNbComponents];
    188      fMassFractionVector = new G4double[maxNbComponents];
     190     fAtomsVector        = new G4int   [fArrayLength];
     191     fMassFractionVector = new G4double[fArrayLength];
    189192  }
    190193
    … …  
    195198     fNumberOfComponents = ++fNumberOfElements;
    196199     element->increaseCountUse();
    197   }
    198   else
    199      G4Exception
     200  } else {
     201    G4cerr << "G4Material::AddElement ERROR for " << fName << " nElement= "
     202           <<  fNumberOfElements << G4endl;
     203    G4Exception
    200204    ("ERROR!!! - Attempt to add more than the declared number of elements.");
    201 
     205  }
    202206  // filled.
    203207  if ( G4int(fNumberOfElements) == maxNbComponents ) {     
    … …  
    226230  // initialization
    227231  if (fNumberOfComponents == 0) {
    228     fMassFractionVector = new G4double[50];
    229     fAtomsVector        = new G4int   [50];
     232    fMassFractionVector = new G4double[fArrayLength];
     233    fAtomsVector        = new G4int   [fArrayLength];
    230234  }
    231235
    … …  
    242246      }
    243247      fNumberOfComponents++; 
    244   }
    245   else
    246      G4Exception
     248  } else {
     249    G4cerr << "G4Material::AddElement ERROR for " << fName << " nElement= "
     250           <<  fNumberOfElements << G4endl;
     251    G4Exception
    247252    ("ERROR!!! - Attempt to add more than the declared number of components.");
    248    
     253  }   
     254
    249255  // filled.
    250256  if (G4int(fNumberOfComponents) == maxNbComponents) {
    … …  
    281287  // initialization
    282288  if (fNumberOfComponents == 0) {
    283     fMassFractionVector = new G4double[50];
    284     fAtomsVector        = new G4int   [50];
     289    fMassFractionVector = new G4double[fArrayLength];
     290    fAtomsVector        = new G4int   [fArrayLength];
     291  }
     292
     293  size_t nelm = material->GetNumberOfElements();
     294
     295  // arrays should be extended
     296  if(nelm > 1) {
     297    G4int nold    = fArrayLength;
     298    fArrayLength += nelm - 1;
     299    G4double* v1 = new G4double[fArrayLength];
     300    G4int* i1    = new G4int[fArrayLength];
     301    for(G4int i=0; i<nold; i++) {
     302      v1[i] = fMassFractionVector[i];
     303      i1[i] = fAtomsVector[i];
     304    }
     305    delete [] fAtomsVector;
     306    delete [] fMassFractionVector;
     307    fMassFractionVector = v1;
     308    fAtomsVector = i1;
    285309  }
    286310
    287311  // filling ...
    288312  if (G4int(fNumberOfComponents) < maxNbComponents) {
    289      for (size_t elm=0; elm < material->GetNumberOfElements(); elm++)
     313     for (size_t elm=0; elm<nelm; elm++)
    290314       {
    291315        G4Element* element = (*(material->GetElementVector()))[elm];
    … …  
    303327       }
    304328      fNumberOfComponents++; 
    305   }
    306   else
    307      G4Exception
    308     ("ERROR!!! - Attempt to add more than the declared number of components.");
    309    
     329  } else {
     330    G4cerr << "G4Material::AddElement ERROR for " << fName << " nElement= "
     331           <<  fNumberOfElements << G4endl;
     332    G4Exception
     333      ("ERROR!!! - Attempt to add more than the declared number of components.");
     334  }   
     335
    310336  // filled.
    311337  if (G4int(fNumberOfComponents) == maxNbComponents) {
    … …  
    541567 
    542568  flux
    543     << " Material: "      << std::setw(8) <<  material->fName
     569    << " Material: "         << std::setw(8) <<  material->fName
    544570    << " " << material->fChemicalFormula << " "
    545     << "  density: "     << std::setw(6) << std::setprecision(3) 
     571    << "  density: "         << std::setw(6) << std::setprecision(3) 
    546572    << G4BestUnit(material->fDensity,"Volumic Mass")
    547     << "  RadL: "        << std::setw(7)  << std::setprecision(3) 
     573    << "  RadL: "            << std::setw(7)  << std::setprecision(3) 
    548574    << G4BestUnit(material->fRadlen,"Length")
    549     << "  Imean: "       << std::setw(7)  << std::setprecision(3) 
     575    << "  Nucl.Int.Length: " << std::setw(7)  << std::setprecision(3) 
     576    << G4BestUnit(material->fNuclInterLen,"Length")   
     577    << "  Imean: "           << std::setw(7)  << std::setprecision(3) 
    550578    << G4BestUnit(material->GetIonisation()->GetMeanExcitationEnergy(),"Energy");
     579   
    551580  if(material->fState == kStateGas)
    552581    flux

  • trunk/source/materials/src/G4MaterialPropertiesTable.cc

    r822 r850  
    2525//
    2626//
    27 // $Id: G4MaterialPropertiesTable.cc,v 1.19 2007/04/25 15:46:55 gum Exp $
    28 // GEANT4 tag $Name: $
     27// $Id: G4MaterialPropertiesTable.cc,v 1.20 2008/06/05 23:38:34 gum Exp $
     28// GEANT4 tag $Name: HEAD $
    2929//
    3030//
    … …  
    8686
    8787void G4MaterialPropertiesTable::AddProperty(const char     *key,
    88                                             G4double *PhotonMomenta,
     88                                            G4double *PhotonEnergies,
    8989                                            G4double *PropertyValues,
    9090                                            G4int     NumEntries)
    … …  
    9494
    9595        G4MaterialPropertyVector *mpv =
    96                         new G4MaterialPropertyVector(PhotonMomenta,
     96                        new G4MaterialPropertyVector(PhotonEnergies,
    9797                                                     PropertyValues,
    9898                                                     NumEntries);
    … …  
    165165
    166166void G4MaterialPropertiesTable::AddEntry(const char     *key,
    167                                          G4double  aPhotonMomentum,
     167                                         G4double  aPhotonEnergy,
    168168                                         G4double  aPropertyValue)
    169169
    … …  
    174174        G4MaterialPropertyVector *targetVector=MPT [G4String(key)];
    175175        if (targetVector != 0) {
    176                 targetVector->AddElement(aPhotonMomentum, aPropertyValue);
     176                targetVector->AddElement(aPhotonEnergy, aPropertyValue);
    177177        }
    178178        else {
    … …  
    183183
    184184void G4MaterialPropertiesTable::RemoveEntry(const char *key, 
    185                                             G4double  aPhotonMomentum)
     185                                            G4double  aPhotonEnergy)
    186186{
    187187//      Allows to remove an entry pair directly from the Material Property Vector
    … …  
    190190        G4MaterialPropertyVector *targetVector=MPT [G4String(key)];
    191191        if (targetVector) {
    192                 targetVector->RemoveElement(aPhotonMomentum);
     192                targetVector->RemoveElement(aPhotonEnergy);
    193193        }
    194194        else {
    … …  
    238238  // fill GROUPVEL vector using RINDEX values
    239239  // rindex built-in "iterator" was advanced to first entry above
    240   G4double E0 = rindex->GetPhotonMomentum();
     240  G4double E0 = rindex->GetPhotonEnergy();
    241241  G4double n0 = rindex->GetProperty();
    242242
    … …  
    247247    // good, we have at least two entries in RINDEX
    248248    // get next energy/value pair
    249     G4double E1 = rindex->GetPhotonMomentum();
     249    G4double E1 = rindex->GetPhotonEnergy();
    250250    G4double n1 = rindex->GetProperty();
    251251
    … …  
    268268      E0 = E1;
    269269      n0 = n1;
    270       E1 = rindex->GetPhotonMomentum();
     270      E1 = rindex->GetPhotonEnergy();
    271271      n1 = rindex->GetProperty();
    272272

  • trunk/source/materials/src/G4MaterialPropertyVector.cc

    r822 r850  
    2525//
    2626//
    27 // $Id: G4MaterialPropertyVector.cc,v 1.15 2006/06/29 19:12:56 gunter Exp $
    28 // GEANT4 tag $Name: $
     27// $Id: G4MaterialPropertyVector.cc,v 1.16 2008/06/05 23:38:51 gum Exp $
     28// GEANT4 tag $Name: HEAD $
    2929//
    3030//
    … …  
    9898        /////////////////
    9999
    100 G4MaterialPropertyVector::G4MaterialPropertyVector(G4double *PhotonMomenta,
     100G4MaterialPropertyVector::G4MaterialPropertyVector(G4double *PhotonEnergies,
    101101                                                   G4double *PropertyValues,
    102102                                                   G4int     NumElements)
    … …  
    106106
    107107        // create a vector filling it with the values
    108         // from PhotonMomenta[] and PropertyValues[]
     108        // from PhotonEnergies[] and PropertyValues[]
    109109
    110110        for(G4int i = 0; i < NumElements; i++) {
    111                 AddElement(PhotonMomenta[i], PropertyValues[i]);
     111                AddElement(PhotonEnergies[i], PropertyValues[i]);
    112112        }
    113113}
    … …  
    141141        // Methods
    142142        ////////////
    143 void G4MaterialPropertyVector::AddElement(G4double aPhotonMomentum,
     143void G4MaterialPropertyVector::AddElement(G4double aPhotonEnergy,
    144144                                          G4double aPropertyValue)
    145145{
    146146        G4MPVEntry *newElement;
    147147       
    148         newElement = new G4MPVEntry(aPhotonMomentum, aPropertyValue);
     148        newElement = new G4MPVEntry(aPhotonEnergy, aPropertyValue);
    149149        MPV.push_back(newElement);
    150150        std::sort(MPV.begin(), MPV.end(), MPVEntry_less());
    … …  
    152152}
    153153
    154 void G4MaterialPropertyVector::RemoveElement(G4double aPhotonMomentum)
     154void G4MaterialPropertyVector::RemoveElement(G4double aPhotonEnergy)
    155155{
    156156        G4MPVEntry *newElement;
    157157        G4MPVEntry *success=0;
    158158
    159         newElement = new G4MPVEntry(aPhotonMomentum, DBL_MAX);
     159        newElement = new G4MPVEntry(aPhotonEnergy, DBL_MAX);
    160160
    161161        std::vector<G4MPVEntry*>::iterator i;
    … …  
    179179
    180180G4double
    181 G4MaterialPropertyVector::GetProperty(G4double aPhotonMomentum) const
     181G4MaterialPropertyVector::GetProperty(G4double aPhotonEnergy) const
    182182{
    183183        G4MPVEntry *target, *temp;
    … …  
    189189        /////////////////////////
    190190
    191         G4double PMmin   = MPV.front()->GetPhotonMomentum();
     191        G4double PMmin   = MPV.front()->GetPhotonEnergy();
    192192        G4double minProp = MPV.front()->GetProperty();
    193         G4double PMmax   = MPV.back() ->GetPhotonMomentum();
     193        G4double PMmax   = MPV.back() ->GetPhotonEnergy();
    194194        G4double maxProp = MPV.back() ->GetProperty();
    195195
    … …  
    198198        ///////////////////////////////////////////
    199199
    200         if (aPhotonMomentum < PMmin)
     200        if (aPhotonEnergy < PMmin)
    201201        {
    202202                G4cout << "\nWarning: G4MaterialPropertyVector::GetProperty";
    … …  
    206206        }
    207207
    208         if (aPhotonMomentum > PMmax)
     208        if (aPhotonEnergy > PMmax)
    209209        {
    210210                G4cout << "\nWarning: G4MaterialPropertyVector::GetProperty";
    … …  
    214214        }
    215215       
    216         target = new G4MPVEntry(aPhotonMomentum, 0.0);
     216        target = new G4MPVEntry(aPhotonEnergy, 0.0);
    217217
    218218        temp = 0;
    … …  
    236236                //////////////////////////////
    237237
    238                 GetAdjacentBins(aPhotonMomentum, &left, &right);
    239 
    240                 pmleft = MPV[left]->GetPhotonMomentum();
    241                 pmright = MPV[right]->GetPhotonMomentum();
    242                 ratio1 = (aPhotonMomentum-pmleft)/(pmright-pmleft);
     238                GetAdjacentBins(aPhotonEnergy, &left, &right);
     239
     240                pmleft = MPV[left]->GetPhotonEnergy();
     241                pmright = MPV[right]->GetPhotonEnergy();
     242                ratio1 = (aPhotonEnergy-pmleft)/(pmright-pmleft);
    243243                ratio2 = 1 - ratio1;
    244244                InterpolatedValue = MPV[left]->GetProperty()*ratio2 +
    … …  
    264264}
    265265
    266 G4double G4MaterialPropertyVector::GetPhotonMomentum() const
     266G4double G4MaterialPropertyVector::GetPhotonEnergy() const
    267267{
    268268//      For use with G4MaterialPropertyVector iterator
    269269
    270270        if(CurrentEntry == -1 || CurrentEntry >= NumEntries) {
    271                 G4Exception("G4MaterialPropertyVector::GetPhotonMomentum ==>"
    272                 "Iterator attempted to Retrieve Photon Momentum out of range");
     271                G4Exception("G4MaterialPropertyVector::GetPhotonEnergy ==>"
     272                "Iterator attempted to Retrieve Photon Energy out of range");
    273273                return DBL_MAX;
    274274        }
    275275        else {
    276                 return MPV[CurrentEntry]->GetPhotonMomentum();
     276                return MPV[CurrentEntry]->GetPhotonEnergy();
    277277        }
    278278}
    279279
    280280G4double
    281 G4MaterialPropertyVector::GetPhotonMomentum(G4double aProperty) const
     281G4MaterialPropertyVector::GetPhotonEnergy(G4double aProperty) const
    282282{
    283283//                              ***NB***
    284 // Assumes that the property is an increasing function of photon momentum (e.g.
     284// Assumes that the property is an increasing function of photon energy (e.g.
    285285// refraction index)
    286286//                              ***NB***
    287287//
    288 // Returns the photon momentum corresponding to the property value passed in.
    289 // If several photon momentum values correspond to the value passed in, the
    290 // function returns the first photon momentum in the vector that corresponds
     288// Returns the photon energy corresponding to the property value passed in.
     289// If several photon energy values correspond to the value passed in, the
     290// function returns the first photon energy in the vector that corresponds
    291291// to that value.
    292292
    … …  
    299299
    300300        G4double PropMin = MPV.front()->GetProperty();
    301         G4double PMmin   = MPV.front()->GetPhotonMomentum();
     301        G4double PMmin   = MPV.front()->GetPhotonEnergy();
    302302        G4double PropMax = MPV.back() ->GetProperty();
    303         G4double PMmax   = MPV.back() ->GetPhotonMomentum();
     303        G4double PMmax   = MPV.back() ->GetPhotonEnergy();
    304304
    305305        ///////////////////////////////////////////
    … …  
    309309        if (aProperty < PropMin)
    310310        {
    311           G4cout << "\nWarning: G4MaterialPropertyVector::GetPhotonMomentum";
    312           G4cout << "\n==> attempt to Retrieve Photon Momentum out of range"
     311          G4cout << "\nWarning: G4MaterialPropertyVector::GetPhotonEnergy";
     312          G4cout << "\n==> attempt to Retrieve Photon Energy out of range"
    313313               << G4endl;
    314314          return PMmin;
    … …  
    316316
    317317        if (aProperty > PropMax) {
    318           G4cout << "\nWarning: G4MaterialPropertyVector::GetPhotonMomentum";
    319           G4cout << "\n==> attempt to Retrieve Photon Momentum out of range"
     318          G4cout << "\nWarning: G4MaterialPropertyVector::GetPhotonEnergy";
     319          G4cout << "\n==> attempt to Retrieve Photon Energy out of range"
    320320               << G4endl;
    321321          return PMmax;
    … …  
    335335                if (MPV[mid]->GetProperty() == aProperty) {
    336336
    337                         // Get first photon momentum value in vector that
     337                        // Get first photon energy value in vector that
    338338                        // corresponds to property value 
    339339
    … …  
    343343                        }
    344344
    345                         InterpolatedValue = MPV[mid]->GetPhotonMomentum();
    346                         goto end_GetPhotonMomentum;     
     345                        InterpolatedValue = MPV[mid]->GetPhotonEnergy();
     346                        goto end_GetPhotonEnergy;       
    347347                }
    348348                if (MPV[mid]->GetProperty() < aProperty)
    … …  
    357357        ratio1 = (aProperty - pleft) / (pright - pleft);
    358358        ratio2 = 1 - ratio1;
    359         InterpolatedValue = MPV[left]->GetPhotonMomentum()*ratio2 +
    360                             MPV[right]->GetPhotonMomentum()*ratio1;
    361 
    362 end_GetPhotonMomentum:
     359        InterpolatedValue = MPV[left]->GetPhotonEnergy()*ratio2 +
     360                            MPV[right]->GetPhotonEnergy()*ratio1;
     361
     362end_GetPhotonEnergy:
    363363
    364364        return InterpolatedValue;
    … …  
    392392
    393393void
    394 G4MaterialPropertyVector::GetAdjacentBins(G4double aPhotonMomentum,
     394G4MaterialPropertyVector::GetAdjacentBins(G4double aPhotonEnergy,
    395395                                          G4int *left,
    396396                                          G4int *right) const
    … …  
    401401        *right = (MPV.size() - 1); // was .entries()
    402402
    403         // find values in bins on either side of aPhotonMomentum
     403        // find values in bins on either side of aPhotonEnergy
    404404
    405405        do {
    406406                mid = (*left + *right)/2;
    407407
    408                 if (MPV[mid]->GetPhotonMomentum() < aPhotonMomentum)
     408                if (MPV[mid]->GetPhotonEnergy() < aPhotonEnergy)
    409409                {
    410410                        *left = mid;

  • trunk/source/materials/src/G4NistElementBuilder.cc

    r822 r850  
    2424// ********************************************************************
    2525//
    26 // $Id: G4NistElementBuilder.cc,v 1.16 2007/07/28 15:58:03 vnivanch Exp $
    27 // GEANT4 tag $Name: $
     26// $Id: G4NistElementBuilder.cc,v 1.22 2008/08/11 11:53:11 vnivanch Exp $
     27// GEANT4 tag $Name: HEAD $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    5858#include "G4NistElementBuilder.hh"
    5959#include "G4Element.hh"
     60#include <sstream>
     61
    6062
    6163//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    … …  
    7779//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    7880
     81G4int G4NistElementBuilder::GetZ(const G4String& name)
     82{
     83  G4int Z = maxNumElements;
     84  do {Z--;} while( Z>0 && elmSymbol[Z] != name);
     85  return Z;
     86}
     87
     88//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
     89
    7990G4Element* G4NistElementBuilder::FindOrBuildElement(const G4String& symb,
    80                                                           G4bool buildIsotopes)
     91                                                    G4bool buildIsotopes)
    8192{
    8293  if(first) {
    … …  
    140151    for (G4int i=0; i<nc; i++) {
    141152       if (relAbundance[idx + i] > 0.0) {
    142          ist = new G4Isotope(elmSymbol[Z],Z, n0 + i,
    143                              massIsotopes[idx + i]*gram/mole);
     153         std::ostringstream os;
     154         os << elmSymbol[Z] << n0 + i;
     155         ist = new G4Isotope(os.str(), Z, n0 + i,
     156                             GetAtomicMass(Z, n0 + i)*g/(mole*amu_c2));
     157         /*
     158         G4cout << " Z= " << Z << " N= " << n0 + i
     159                << " miso(amu)= " <<  GetIsotopeMass(Z, n0 + i)/amu_c2
     160                << " matom(amu)= " << GetAtomicMass(Z, n0 + i)/amu_c2 << G4endl;
     161         */
    144162         iso.push_back(ist);
    145163       }
    … …  
    193211      G4cout << G4endl;
    194212      G4cout << "          mass(amu): ";
    195       for(j=0; j<nc; j++) {G4cout << massIsotopes[idx + j] << " ";}
     213      for(j=0; j<nc; j++) {G4cout << GetAtomicMass(i, n0 + j) << " ";}
    196214      G4cout << G4endl;
    197215      G4cout << "     abanbance: ";
    … …  
    232250  G4double www;
    233251  size_t nm = nc;
    234   //  G4double delm = G4double(Z)*electron_mass_c2/amu_c2;
    235252
    236253  for(size_t i=0; i<nm; i++) {
    237254    www = 0.01*(&W)[i];
    238     massIsotopes[index] = (&A)[i]; // - delm;
    239     sigMass[index]      = (&sA)[i];
     255    // mass of the isotope in G4 units
     256    massIsotopes[index] = (&A)[i]*amu_c2 - Z*electron_mass_c2 + bindingEnergy[Z];
     257    sigMass[index]      = (&sA)[i]*amu_c2;
    240258    relAbundance[index] = www;
    241259
     260    // computation of mean atomic mass of the element in atomic units
    242261    atomicMass[Z] += www*(&A)[i];
    243262    ww += www;
    … …  
    258277void G4NistElementBuilder::Initialise()
    259278{
     279  // Parameterisation from D.Lunney,J.M.Pearson,C.Thibault,
     280  // Rev.Mod.Phys. 75 (2003) 1021
     281  bindingEnergy[0] = 0.0;
     282  for(G4int i=1; i<maxNumElements; i++) {
     283    G4double Z = G4double(i);
     284    bindingEnergy[i] = (14.4381*std::pow(Z,2.39) + 1.55468e-6*std::pow(Z,5.35))*eV;
     285  }
     286
     287  // NIST data
    260288  index    = 0;
    261289   
    … …  
    266294  double HA[6] =
    267295  {1.00783, 2.0141, 3.01605, 4.02783, 5.03954, 6.04494};
     296
     297  // Garantee consistence with G4 masses
     298  HA[0] = (proton_mass_c2 + electron_mass_c2 - bindingEnergy[1])/amu_c2;
     299  HA[1] = (1.875613*GeV   + electron_mass_c2 - bindingEnergy[1])/amu_c2;
     300  HA[2] = (2.80925*GeV    + electron_mass_c2 - bindingEnergy[1])/amu_c2;
    268301 
    269302  double HS[6] =
    … …  
    281314  double HeA[8] =
    282315  {3.01603, 4.0026, 5.01222, 6.01889, 7.02803, 8.03392, 9.04382, 10.0524};
     316
     317  // Garantee consistence with G4 masses
     318  HeA[0] = (2.80923*GeV  + 2.0*electron_mass_c2 - bindingEnergy[2])/amu_c2;
     319  HeA[1] = (3.727417*GeV + 2.0*electron_mass_c2 - bindingEnergy[2])/amu_c2;
    283320 
    284321  double HeS[8] =

  • trunk/source/materials/src/G4NistManager.cc

    r822 r850  
    2424// ********************************************************************
    2525//
    26 // $Id: G4NistManager.cc,v 1.15 2007/12/11 13:32:08 gcosmo Exp $
    27 // GEANT4 tag $Name: geant4-09-01-patch-02 $
     26// $Id: G4NistManager.cc,v 1.19 2008/07/23 14:49:31 vnivanch Exp $
     27// GEANT4 tag $Name: HEAD $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    6363
    6464G4NistManager* G4NistManager::instance = 0;
    65 G4double G4NistManager::POWERZ13[256] = {0};
    66 G4double G4NistManager::LOGA[256] = {0};
    6765
    6866//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
    … …  
    8987 
    9088  messenger  = new G4NistMessenger(this); 
     89
     90  // compute frequently used values
    9191  for(G4int i=1; i<256; i++) {
    9292    G4double x = G4double(i);
    … …  
    9494    LOGA[i] = std::log(x);
    9595  }
     96  for(G4int j=1; j<101; j++) {
     97    G4double A = elmBuilder->GetA(j);
     98    POWERA27[j] = std::pow(A,0.27);
     99    LOGAZ[j]    = std::log(A);
     100  }
    96101  POWERZ13[0] = 1.0;
     102  POWERA27[0] = 1.0;
    97103  LOGA[0]     = 0.0;
     104  LOGAZ[0]    = 0.0;
    98105}
    99106
    … …  
    176183
    177184//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    178 
    179 G4double G4NistManager::GetZ13(G4double Z)
    180 {
    181   G4int iz = G4int(Z);
    182   G4double x = (Z - G4double(iz))/(3.0*Z);
    183   if(iz > 255) iz = 255;
    184   else if(iz < 0) iz = 0;
    185   return POWERZ13[iz]*(1.0 + x);
    186 }
    187 
    188 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    189 
    190 G4double G4NistManager::GetLOGA(G4double A)
    191 {
    192   G4int ia = G4int(A);
    193   G4double x = (A - G4double(ia))/A;
    194   if(ia > 255) ia = 255;
    195   else if(ia < 0) ia = 0;
    196   return LOGA[ia] + x;
    197 }
    198 
    199 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

  • trunk/source/materials/src/G4NistMaterialBuilder.cc

    r822 r850  
    2424// ********************************************************************
    2525//
    26 // $Id: G4NistMaterialBuilder.cc,v 1.17.2.3 2008/04/28 08:39:03 gcosmo Exp $
    27 // GEANT4 tag $Name: geant4-09-01-patch-02 $
     26// $Id: G4NistMaterialBuilder.cc,v 1.19 2008/04/28 08:51:29 vnivanch Exp $
     27// GEANT4 tag $Name: HEAD $
    2828//
    2929//
    … …  
    212212    return 0;
    213213  }
     214
     215  // add parameters of material into internal vectors
     216  // density in g/cm3, mean ionisation potential is not defined
    214217  AddMaterial(name,dens*cm3/g,0,0.,nm,state,temp,pressure);
     218
    215219  for (G4int i=0; i<nm; i++) {
    216220    G4int Z = G4int((elmBuilder->FindOrBuildElement(elm[i]))->GetZ());
    … …  
    240244    return 0;
    241245  }
     246
     247  // add parameters of material into internal vectors
     248  // density in g/cm3, mean ionisation potential is not defined
    242249  AddMaterial(name,dens*cm3/g,0,0.,nm,state,temp,pressure);
     250
    243251  for (G4int i=0; i<nm; i++) {
    244252    G4int Z = G4int((elmBuilder->FindOrBuildElement(elm[i]))->GetZ());
    … …  
    418426                                        G4double temp, G4double pres)
    419427{
     428  // add parameters of material into internal vectors
     429  // density in g/cm3, mean ionisation potential in eV
     430
    420431  if (nCurrent != 0) {
    421432    G4cout << "WARNING: G4NistMaterialBuilder::AddMaterial problem: previous "

  • trunk/source/materials/src/G4NistMessenger.cc

    r822 r850  
    2626//
    2727// $Id: G4NistMessenger.cc,v 1.4 2007/05/02 10:48:52 vnivanch Exp $
    28 // GEANT4 tag $Name: $
     28// GEANT4 tag $Name: HEAD $
    2929//
    3030//

  • trunk/source/materials/src/G4OpticalSurface.cc

    r822 r850  
    2525//
    2626//
    27 // $Id: G4OpticalSurface.cc,v 1.10 2006/06/29 19:13:08 gunter Exp $
    28 // GEANT4 tag $Name: $
     27// $Id: G4OpticalSurface.cc,v 1.11 2008/07/21 20:55:51 gum Exp $
     28// GEANT4 tag $Name: HEAD $
    2929//
    3030//
    … …  
    104104}
    105105
    106 G4OpticalSurface::~G4OpticalSurface(){}
     106G4OpticalSurface::~G4OpticalSurface()
     107{
     108  Overwrite();
     109}
    107110
    108111G4int G4OpticalSurface::operator==(const G4OpticalSurface &right) const

  • trunk/source/materials/src/G4SandiaTable.cc

    r822 r850  
    2626//
    2727// $Id: G4SandiaTable.cc,v 1.34 2007/10/02 10:13:33 vnivanch Exp $
    28 // GEANT4 tag $Name: $
     28// GEANT4 tag $Name: HEAD $
    2929//
    3030//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....

  • trunk/source/materials/src/G4SurfaceProperty.cc

    r822 r850  
    2525//
    2626//
    27 // $Id: G4SurfaceProperty.cc,v 1.1 2007/04/25 16:19:26 gum Exp $
    28 // GEANT4 tag $Name: $
     27// $Id: G4SurfaceProperty.cc,v 1.4 2008/07/21 20:55:34 gum Exp $
     28// GEANT4 tag $Name: HEAD $
    2929//
    3030//
    … …  
    5252
    5353//
     54// Constructor and destructor
     55//
     56G4SurfaceProperty::G4SurfaceProperty( const G4String& name,
     57                                            G4SurfaceType type )
     58  : theName(name), theType(type)
     59{
     60  theSurfacePropertyTable.push_back(this);
     61}
     62
     63G4SurfaceProperty::~G4SurfaceProperty()
     64{
     65  Overwrite();
     66}
     67
     68//
    5469// Methods
    5570//
    … …  
    6782// Dump info for known surface properties
    6883//
    69 void G4SurfaceProperty::DumpInfo()
     84void G4SurfaceProperty::DumpTableInfo()
    7085{
    7186  G4cout << "***** Surface Property Table : Nb of Surface Properties = "
    … …  
    85100void G4SurfaceProperty::CleanSurfacePropertyTable()
    86101{
    87   DumpInfo();
     102  DumpTableInfo();
    88103  G4SurfacePropertyTable::iterator pos;
    89104  for(pos=theSurfacePropertyTable.begin();
    … …  
    93108  }
    94109  theSurfacePropertyTable.clear();
    95   DumpInfo();
     110  DumpTableInfo();
    96111}
Note: See TracChangeset for help on using the changeset viewer.