Geant4 Cross Reference

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

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

Diff markup

Differences between /materials/src/G4IonisParamMat.cc (Version 11.3.0) and /materials/src/G4IonisParamMat.cc (Version 10.6.p3)


  1 //                                                  1 //
  2 // *******************************************      2 // ********************************************************************
  3 // * License and Disclaimer                         3 // * License and Disclaimer                                           *
  4 // *                                                4 // *                                                                  *
  5 // * The  Geant4 software  is  copyright of th      5 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
  6 // * the Geant4 Collaboration.  It is provided      6 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
  7 // * conditions of the Geant4 Software License      7 // * conditions of the Geant4 Software License,  included in the file *
  8 // * LICENSE and available at  http://cern.ch/      8 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
  9 // * include a list of copyright holders.           9 // * include a list of copyright holders.                             *
 10 // *                                               10 // *                                                                  *
 11 // * Neither the authors of this software syst     11 // * Neither the authors of this software system, nor their employing *
 12 // * institutes,nor the agencies providing fin     12 // * institutes,nor the agencies providing financial support for this *
 13 // * work  make  any representation or  warran     13 // * work  make  any representation or  warranty, express or implied, *
 14 // * regarding  this  software system or assum     14 // * regarding  this  software system or assume any liability for its *
 15 // * use.  Please see the license in the file      15 // * use.  Please see the license in the file  LICENSE  and URL above *
 16 // * for the full disclaimer and the limitatio     16 // * for the full disclaimer and the limitation of liability.         *
 17 // *                                               17 // *                                                                  *
 18 // * This  code  implementation is the result      18 // * This  code  implementation is the result of  the  scientific and *
 19 // * technical work of the GEANT4 collaboratio     19 // * technical work of the GEANT4 collaboration.                      *
 20 // * By using,  copying,  modifying or  distri     20 // * By using,  copying,  modifying or  distributing the software (or *
 21 // * any work based  on the software)  you  ag     21 // * any work based  on the software)  you  agree  to acknowledge its *
 22 // * use  in  resulting  scientific  publicati     22 // * use  in  resulting  scientific  publications,  and indicate your *
 23 // * acceptance of all terms of the Geant4 Sof     23 // * acceptance of all terms of the Geant4 Software license.          *
 24 // *******************************************     24 // ********************************************************************
                                                   >>  25 //
                                                   >>  26 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
 25                                                    27 
 26 // 09-07-98, data moved from G4Material, M.Mai     28 // 09-07-98, data moved from G4Material, M.Maire
 27 // 18-07-98, bug corrected in ComputeDensityEf     29 // 18-07-98, bug corrected in ComputeDensityEffect() for gas
 28 // 16-01-01, bug corrected in ComputeDensityEf     30 // 16-01-01, bug corrected in ComputeDensityEffect() E100eV (L.Urban)
 29 // 08-02-01, fShellCorrectionVector correctly      31 // 08-02-01, fShellCorrectionVector correctly handled (mma)
 30 // 28-10-02, add setMeanExcitationEnergy (V.Iv     32 // 28-10-02, add setMeanExcitationEnergy (V.Ivanchenko)
 31 // 06-09-04, factor 2 to shell correction term <<  33 // 06-09-04, factor 2 to shell correction term (V.Ivanchenko) 
 32 // 10-05-05, add a missing coma in FindMeanExc     34 // 10-05-05, add a missing coma in FindMeanExcitationEnergy() - Bug#746 (mma)
 33 // 27-09-07, add computation of parameters for     35 // 27-09-07, add computation of parameters for ions (V.Ivanchenko)
 34 // 04-03-08, remove reference to G4NistManager     36 // 04-03-08, remove reference to G4NistManager. Add fBirks constant (mma)
 35 // 30-10-09, add G4DensityEffectData class and     37 // 30-10-09, add G4DensityEffectData class and density effect computation (VI)
 36                                                    38 
 37 #include "G4IonisParamMat.hh"                  <<  39 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
 38                                                    40 
 39 #include "G4AtomicShells.hh"                   <<  41 #include "G4IonisParamMat.hh"
 40 #include "G4AutoLock.hh"                       << 
 41 #include "G4DensityEffectData.hh"              << 
 42 #include "G4Exp.hh"                            << 
 43 #include "G4Log.hh"                            << 
 44 #include "G4Material.hh"                           42 #include "G4Material.hh"
                                                   >>  43 #include "G4DensityEffectData.hh"
                                                   >>  44 #include "G4DensityEffectCalculator.hh"
                                                   >>  45 #include "G4AtomicShells.hh"
 45 #include "G4NistManager.hh"                        46 #include "G4NistManager.hh"
 46 #include "G4PhysicalConstants.hh"              << 
 47 #include "G4Pow.hh"                                47 #include "G4Pow.hh"
                                                   >>  48 #include "G4PhysicalConstants.hh"
 48 #include "G4SystemOfUnits.hh"                      49 #include "G4SystemOfUnits.hh"
 49                                                    50 
 50 G4DensityEffectData* G4IonisParamMat::fDensity     51 G4DensityEffectData* G4IonisParamMat::fDensityData = nullptr;
 51                                                    52 
 52 namespace                                      <<  53 #ifdef G4MULTITHREADED
 53 {                                              <<  54   G4Mutex G4IonisParamMat::ionisMutex = G4MUTEX_INITIALIZER;
 54   G4Mutex ionisMutex = G4MUTEX_INITIALIZER;    <<  55 #endif
 55 }                                              << 
 56                                                    56 
 57 //....oooOO0OOooo........oooOO0OOooo........oo     57 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
 58                                                    58 
 59 G4IonisParamMat::G4IonisParamMat(const G4Mater <<  59 G4IonisParamMat::G4IonisParamMat(const G4Material* material)
                                                   >>  60   : fMaterial(material)
 60 {                                                  61 {
 61   fBirks = 0.;                                     62   fBirks = 0.;
 62   fMeanEnergyPerIon = 0.0;                         63   fMeanEnergyPerIon = 0.0;
 63   twoln10 = 2. * G4Pow::GetInstance()->logZ(10 <<  64   twoln10 = 2.*G4Pow::GetInstance()->logZ(10);
 64                                                    65 
 65   // minimal set of default parameters for den     66   // minimal set of default parameters for density effect
 66   fCdensity = 0.0;                                 67   fCdensity = 0.0;
 67   fD0density = 0.0;                                68   fD0density = 0.0;
 68   fAdjustmentFactor = 1.0;                         69   fAdjustmentFactor = 1.0;
 69   if (fDensityData == nullptr) {               <<  70   if(fDensityData == nullptr) { fDensityData = new G4DensityEffectData(); }
 70     fDensityData = new G4DensityEffectData();  << 
 71   }                                            << 
 72   fDensityEffectCalc = nullptr;                    71   fDensityEffectCalc = nullptr;
 73                                                    72 
 74   // compute parameters                            73   // compute parameters
 75   ComputeMeanParameters();                         74   ComputeMeanParameters();
 76   ComputeDensityEffectParameters();                75   ComputeDensityEffectParameters();
 77   ComputeFluctModel();                             76   ComputeFluctModel();
 78   ComputeIonParameters();                          77   ComputeIonParameters();
 79 }                                                  78 }
 80                                                    79 
 81 //....oooOO0OOooo........oooOO0OOooo........oo     80 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
 82                                                    81 
 83 G4IonisParamMat::~G4IonisParamMat()            <<  82 // Fake default constructor - sets only member data and allocates memory
                                                   >>  83 //                            for usage restricted to object persistency
                                                   >>  84 
                                                   >>  85 G4IonisParamMat::G4IonisParamMat(__void__&)
                                                   >>  86   : fMaterial(nullptr), fShellCorrectionVector(nullptr)
 84 {                                                  87 {
 85   delete fDensityEffectCalc;                   <<  88   fMeanExcitationEnergy = 0.0;
 86   delete[] fShellCorrectionVector;             <<  89   fLogMeanExcEnergy = 0.0;
 87   delete fDensityData;                         <<  90   fTaul = 0.0;
 88   fDensityData = nullptr;                      <<  91   fCdensity = 0.0;
 89   fShellCorrectionVector = nullptr;            <<  92   fMdensity = 0.0;
                                                   >>  93   fAdensity = 0.0;
                                                   >>  94   fX0density = 0.0;
                                                   >>  95   fX1density = 0.0;
                                                   >>  96   fD0density = 0.0;
                                                   >>  97   fPlasmaEnergy = 0.0;
                                                   >>  98   fAdjustmentFactor = 0.0;
                                                   >>  99   fF1fluct = 0.0;          
                                                   >> 100   fF2fluct = 0.0;                       
                                                   >> 101   fEnergy1fluct = 0.0;
                                                   >> 102   fLogEnergy1fluct = 0.0;
                                                   >> 103   fEnergy2fluct = 0.0;
                                                   >> 104   fLogEnergy2fluct = 0.0;
                                                   >> 105   fEnergy0fluct = 0.0;
                                                   >> 106   fRateionexcfluct = 0.0;
                                                   >> 107   fZeff = 0.0;
                                                   >> 108   fFermiEnergy = 0.0;
                                                   >> 109   fLfactor = 0.0;
                                                   >> 110   fInvA23 = 0.0;
                                                   >> 111   fBirks = 0.0;
                                                   >> 112   fMeanEnergyPerIon = 0.0;
                                                   >> 113   twoln10 = 2.*G4Pow::GetInstance()->logZ(10);
                                                   >> 114 
 90   fDensityEffectCalc = nullptr;                   115   fDensityEffectCalc = nullptr;
                                                   >> 116   if(fDensityData == nullptr) { fDensityData = new G4DensityEffectData(); }
 91 }                                                 117 }
 92                                                   118 
 93 //....oooOO0OOooo........oooOO0OOooo........oo    119 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
 94                                                   120 
 95 G4double G4IonisParamMat::GetDensityCorrection << 121 G4IonisParamMat::~G4IonisParamMat()
 96 {                                                 122 {
 97   // x = log10(beta*gamma)                     << 123   delete fDensityEffectCalc;
 98   G4double y = 0.0;                            << 124   delete [] fShellCorrectionVector; 
 99   if (x < fX0density) {                        << 125   delete fDensityData;
100     if (fD0density > 0.0) {                    << 126   fDensityData = nullptr;
101       y = fD0density * G4Exp(twoln10 * (x - fX << 127   fShellCorrectionVector = nullptr;
102     }                                          << 128   fDensityEffectCalc = nullptr;
103   }                                            << 
104   else if (x >= fX1density) {                  << 
105     y = twoln10 * x - fCdensity;               << 
106   }                                            << 
107   else {                                       << 
108     y = twoln10 * x - fCdensity + fAdensity *  << 
109   }                                            << 
110   return y;                                    << 
111 }                                                 129 }
112                                                   130 
113 //....oooOO0OOooo........oooOO0OOooo........oo    131 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
114                                                   132 
115 void G4IonisParamMat::ComputeMeanParameters()     133 void G4IonisParamMat::ComputeMeanParameters()
116 {                                                 134 {
117   // compute mean excitation energy and shell     135   // compute mean excitation energy and shell correction vector
118   fTaul = (*(fMaterial->GetElementVector()))[0    136   fTaul = (*(fMaterial->GetElementVector()))[0]->GetIonisation()->GetTaul();
119                                                   137 
120   std::size_t nElements = fMaterial->GetNumber << 138   fMeanExcitationEnergy = 0.;
                                                   >> 139   fLogMeanExcEnergy = 0.;
                                                   >> 140 
                                                   >> 141   size_t nElements = fMaterial->GetNumberOfElements();
121   const G4ElementVector* elmVector = fMaterial    142   const G4ElementVector* elmVector = fMaterial->GetElementVector();
122   const G4double* nAtomsPerVolume = fMaterial-    143   const G4double* nAtomsPerVolume = fMaterial->GetVecNbOfAtomsPerVolume();
123                                                << 144  
124   fMeanExcitationEnergy = FindMeanExcitationEn    145   fMeanExcitationEnergy = FindMeanExcitationEnergy(fMaterial);
125   fLogMeanExcEnergy = 0.;                      << 
126                                                   146 
127   // Chemical formula defines mean excitation     147   // Chemical formula defines mean excitation energy
128   if (fMeanExcitationEnergy > 0.0) {           << 148   if(fMeanExcitationEnergy > 0.0) {
129     fLogMeanExcEnergy = G4Log(fMeanExcitationE    149     fLogMeanExcEnergy = G4Log(fMeanExcitationEnergy);
130                                                   150 
131     // Compute average                         << 151     // Compute average 
132   }                                            << 152   } else {
133   else {                                       << 153     for (size_t i=0; i < nElements; i++) {
134     for (std::size_t i = 0; i < nElements; ++i << 
135       const G4Element* elm = (*elmVector)[i];     154       const G4Element* elm = (*elmVector)[i];
136       fLogMeanExcEnergy +=                     << 155       fLogMeanExcEnergy += nAtomsPerVolume[i]*elm->GetZ()
137         nAtomsPerVolume[i] * elm->GetZ() * G4L << 156   *G4Log(elm->GetIonisation()->GetMeanExcitationEnergy());
138     }                                             157     }
139     fLogMeanExcEnergy /= fMaterial->GetTotNbOf    158     fLogMeanExcEnergy /= fMaterial->GetTotNbOfElectPerVolume();
140     fMeanExcitationEnergy = G4Exp(fLogMeanExcE    159     fMeanExcitationEnergy = G4Exp(fLogMeanExcEnergy);
141   }                                               160   }
142                                                   161 
143   fShellCorrectionVector = new G4double[3];    << 162   fShellCorrectionVector = new G4double[3]; 
144                                                   163 
145   for (G4int j = 0; j <= 2; ++j) {             << 164   for (G4int j=0; j<=2; j++)
                                                   >> 165   {
146     fShellCorrectionVector[j] = 0.;               166     fShellCorrectionVector[j] = 0.;
147                                                   167 
148     for (std::size_t k = 0; k < nElements; ++k << 168     for (size_t k=0; k<nElements; k++) {
149       fShellCorrectionVector[j] +=             << 169       fShellCorrectionVector[j] += nAtomsPerVolume[k]
150         nAtomsPerVolume[k] * (((*elmVector)[k] << 170   *(((*elmVector)[k])->GetIonisation()->GetShellCorrectionVector())[j];
151     }                                             171     }
152     fShellCorrectionVector[j] *= 2.0 / fMateri << 172     fShellCorrectionVector[j] *= 2.0/fMaterial->GetTotNbOfElectPerVolume();
153   }                                            << 173   } 
154 }                                                 174 }
155                                                   175 
156 //....oooOO0OOooo........oooOO0OOooo........oo    176 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
157                                                   177 
158 G4DensityEffectData* G4IonisParamMat::GetDensi << 178 G4DensityEffectData* G4IonisParamMat::GetDensityEffectData()
                                                   >> 179 {
159   return fDensityData;                            180   return fDensityData;
160 }                                                 181 }
161                                                   182 
162 //....oooOO0OOooo........oooOO0OOooo........oo    183 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
163                                                   184 
                                                   >> 185 G4double G4IonisParamMat::DensityCorrection(G4double x)
                                                   >> 186 {
                                                   >> 187   return (!fDensityEffectCalc) ? GetDensityCorrection(x) 
                                                   >> 188     : fDensityEffectCalc->ComputeDensityCorrection(x);
                                                   >> 189 }
                                                   >> 190 
                                                   >> 191 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
                                                   >> 192                     
164 void G4IonisParamMat::ComputeDensityEffectPara    193 void G4IonisParamMat::ComputeDensityEffectParameters()
165 {                                                 194 {
166   G4State State = fMaterial->GetState();          195   G4State State = fMaterial->GetState();
167   G4double density = fMaterial->GetDensity();  << 
168                                                   196 
169   // Check if density effect data exist in the    197   // Check if density effect data exist in the table
170   // R.M. Sternheimer, Atomic Data and Nuclear    198   // R.M. Sternheimer, Atomic Data and Nuclear Data Tables, 30: 261 (1984)
171   // or is assign to one of data set in this t    199   // or is assign to one of data set in this table
172   G4int idx = fDensityData->GetIndex(fMaterial    200   G4int idx = fDensityData->GetIndex(fMaterial->GetName());
173   auto nelm = (G4int)fMaterial->GetNumberOfEle << 201   G4int nelm= fMaterial->GetNumberOfElements();
174   G4int Z0 = ((*(fMaterial->GetElementVector() << 202   G4int Z0  = ((*(fMaterial->GetElementVector()))[0])->GetZasInt();
175   const G4Material* bmat = fMaterial->GetBaseM    203   const G4Material* bmat = fMaterial->GetBaseMaterial();
176   G4NistManager* nist = G4NistManager::Instanc    204   G4NistManager* nist = G4NistManager::Instance();
177                                                   205 
178   // arbitrary empirical limits                << 206   // arbitrary empirical limits 
179   // parameterisation with very different dens    207   // parameterisation with very different density is not applicable
180   static const G4double corrmax = 1.;             208   static const G4double corrmax = 1.;
181   static const G4double massfracmax = 0.9;        209   static const G4double massfracmax = 0.9;
182                                                   210 
183   // for simple non-NIST materials             << 211   // for simple non-NIST materials 
184   G4double corr = 0.0;                            212   G4double corr = 0.0;
185                                                   213 
186   if (idx < 0 && 1 == nelm) {                  << 214   if(idx < 0 && 1 == nelm) {
187     G4int z = (1 == Z0 && State == kStateLiqui << 215     idx = fDensityData->GetElementIndex(Z0, fMaterial->GetState());
188     idx = fDensityData->GetElementIndex(z);    << 
189                                                   216 
190     // Correction for base material or for non    217     // Correction for base material or for non-nominal density
191     // Except cases of very different density     218     // Except cases of very different density defined in user code
192     if (idx >= 0 && 0 < z) {                   << 219     if(idx >= 0) {
193       G4double dens = nist->GetNominalDensity(    220       G4double dens = nist->GetNominalDensity(Z0);
194       if (dens <= 0.0) {                       << 221       if(dens <= 0.0) { idx = -1; }
195         idx = -1;                              << 
196       }                                        << 
197       else {                                      222       else {
198         corr = G4Log(dens / density);          << 223   corr = G4Log(dens/fMaterial->GetDensity());
199         if (std::abs(corr) > corrmax) {        << 
200           idx = -1;                            << 
201         }                                      << 
202       }                                           224       }
                                                   >> 225       if(std::abs(corr) > corrmax) { idx = -1; }
203     }                                             226     }
204   }                                               227   }
205   // for base material case                       228   // for base material case
206   if (idx < 0 && nullptr != bmat) {            << 229   if(idx < 0 && bmat) {
207     idx = fDensityData->GetIndex(bmat->GetName    230     idx = fDensityData->GetIndex(bmat->GetName());
208     if (idx >= 0) {                            << 231     if(idx >= 0) {
209       corr = G4Log(bmat->GetDensity() / densit << 232       corr = G4Log(bmat->GetDensity()/fMaterial->GetDensity());
210       if (std::abs(corr) > corrmax) {          << 233       if(std::abs(corr) > corrmax) { idx = -1; }
211         idx = -1;                              << 
212       }                                        << 
213     }                                             234     }
214   }                                               235   }
215                                                   236 
216   // for compound non-NIST materials with one  << 237   // for compound non-NIST materials with one element dominating 
217   if (idx < 0 && 1 < nelm) {                   << 238   if(idx < 0 && 1 < nelm) {
218     const G4double tot = fMaterial->GetTotNbOf    239     const G4double tot = fMaterial->GetTotNbOfAtomsPerVolume();
219     for (G4int i = 0; i < nelm; ++i) {         << 240     for(G4int i=0; i<nelm; ++i) {
220       const G4double frac = fMaterial->GetVecN << 241       const G4double frac = fMaterial->GetVecNbOfAtomsPerVolume()[i]/tot;
221       if (frac > massfracmax) {                << 242       if(frac > massfracmax) {
222         Z0 = ((*(fMaterial->GetElementVector()    243         Z0 = ((*(fMaterial->GetElementVector()))[i])->GetZasInt();
223         idx = fDensityData->GetElementIndex(Z0 << 244   idx = fDensityData->GetElementIndex(Z0, fMaterial->GetState());
224         G4double dens = nist->GetNominalDensit << 245   G4double dens = nist->GetNominalDensity(Z0);
225         if (idx >= 0 && dens > 0.0) {          << 246         if(idx >= 0 && dens > 0.0) {
226           corr = G4Log(dens / density);        << 247     corr = G4Log(dens/fMaterial->GetDensity());
227           if (std::abs(corr) > corrmax) {      << 248     if(std::abs(corr) > corrmax) { idx = -1; }
228             idx = -1;                          << 249           else { break; }
229           }                                    << 250   }
230           else {                               << 
231             break;                             << 
232           }                                    << 
233         }                                      << 
234       }                                           251       }
235     }                                             252     }
236   }                                               253   }
                                                   >> 254   //G4cout<<"DensityEffect for "<<fMaterial->GetName()<<"  "<< idx << G4endl; 
                                                   >> 255 
                                                   >> 256   if(idx >= 0) {
237                                                   257 
238   if (idx >= 0) {                              << 
239     // Take parameters for the density effect     258     // Take parameters for the density effect correction from
240     // R.M. Sternheimer et al. Density Effect  << 259     // R.M. Sternheimer et al. Density Effect For The Ionization Loss 
241     // of Charged Particles in Various Substan << 260     // of Charged Particles in Various Substances. 
242     // Atom. Data Nucl. Data Tabl. 30 (1984) 2 << 261     // Atom. Data Nucl. Data Tabl. 30 (1984) 261-271. 
243                                                << 262 
244     fCdensity = fDensityData->GetCdensity(idx) << 263     fCdensity  = fDensityData->GetCdensity(idx); 
245     fMdensity = fDensityData->GetMdensity(idx) << 264     fMdensity  = fDensityData->GetMdensity(idx);
246     fAdensity = fDensityData->GetAdensity(idx) << 265     fAdensity  = fDensityData->GetAdensity(idx);
247     fX0density = fDensityData->GetX0density(id    266     fX0density = fDensityData->GetX0density(idx);
248     fX1density = fDensityData->GetX1density(id    267     fX1density = fDensityData->GetX1density(idx);
249     fD0density = fDensityData->GetDelta0densit    268     fD0density = fDensityData->GetDelta0density(idx);
250     fPlasmaEnergy = fDensityData->GetPlasmaEne    269     fPlasmaEnergy = fDensityData->GetPlasmaEnergy(idx);
251     fAdjustmentFactor = fDensityData->GetAdjus    270     fAdjustmentFactor = fDensityData->GetAdjustmentFactor(idx);
252                                                   271 
253     // parameter C is computed and not taken f    272     // parameter C is computed and not taken from Sternheimer tables
254     // fCdensity = 1. + 2*G4Log(fMeanExcitatio << 273     //fCdensity = 1. + 2*G4Log(fMeanExcitationEnergy/fPlasmaEnergy);
255     // G4cout << "IonisParamMat: " << fMateria << 274     //G4cout << "IonisParamMat: " << fMaterial->GetName() 
256     //     << "  Cst= " << Cdensity << " C= "     275     //     << "  Cst= " << Cdensity << " C= " << fCdensity << G4endl;
257                                                   276 
258     // correction on nominal density              277     // correction on nominal density
259     fCdensity += corr;                         << 278     fCdensity  += corr;
260     fX0density += corr / twoln10;              << 279     fX0density += corr/twoln10;
261     fX1density += corr / twoln10;              << 280     fX1density += corr/twoln10;
262   }                                            << 281 
263   else {                                       << 282   } else {
264     static const G4double Cd2 = 4 * CLHEP::pi  << 283 
265     fPlasmaEnergy = std::sqrt(Cd2 * fMaterial- << 284     static const G4double Cd2 = 4*pi*hbarc_squared*classic_electr_radius;
                                                   >> 285     fPlasmaEnergy = std::sqrt(Cd2*fMaterial->GetTotNbOfElectPerVolume());
266                                                   286 
267     // Compute parameters for the density effe    287     // Compute parameters for the density effect correction in DE/Dx formula.
268     // The parametrization is from R.M. Sternh    288     // The parametrization is from R.M. Sternheimer, Phys. Rev.B,3:3681 (1971)
269     G4int icase;                                  289     G4int icase;
270                                                << 290     
271     fCdensity = 1. + 2 * G4Log(fMeanExcitation << 291     fCdensity = 1. + 2*G4Log(fMeanExcitationEnergy/fPlasmaEnergy);
272     //                                            292     //
273     // condensed materials                        293     // condensed materials
274     //                                         << 294     //  
275     if ((State == kStateSolid) || (State == kS << 295     if ((State == kStateSolid)||(State == kStateLiquid)) {
276       static const G4double E100eV = 100. * CL << 
277       static const G4double ClimiS[] = {3.681, << 
278       static const G4double X0valS[] = {1.0, 1 << 
279       static const G4double X1valS[] = {2.0, 3 << 
280                                                << 
281       if (fMeanExcitationEnergy < E100eV) {    << 
282         icase = 0;                             << 
283       }                                        << 
284       else {                                   << 
285         icase = 1;                             << 
286       }                                        << 
287                                                << 
288       if (fCdensity < ClimiS[icase]) {         << 
289         fX0density = 0.2;                      << 
290       }                                        << 
291       else {                                   << 
292         fX0density = 0.326 * fCdensity - X0val << 
293       }                                        << 
294                                                   296 
295       fX1density = X1valS[icase];              << 297       static const G4double E100eV  = 100.*eV; 
296       fMdensity = 3.0;                         << 298       static const G4double ClimiS[] = {3.681 , 5.215 };
297                                                << 299       static const G4double X0valS[] = {1.0   , 1.5   };
298       // special: Hydrogen                     << 300       static const G4double X1valS[] = {2.0   , 3.0   };
                                                   >> 301                                 
                                                   >> 302       if(fMeanExcitationEnergy < E100eV) { icase = 0; }
                                                   >> 303       else                               { icase = 1; } 
                                                   >> 304 
                                                   >> 305       if(fCdensity < ClimiS[icase])    { fX0density = 0.2; }
                                                   >> 306       else { fX0density = 0.326*fCdensity - X0valS[icase]; }
                                                   >> 307 
                                                   >> 308       fX1density = X1valS[icase]; fMdensity = 3.0;
                                                   >> 309       
                                                   >> 310       //special: Hydrogen
299       if (1 == nelm && 1 == Z0) {                 311       if (1 == nelm && 1 == Z0) {
300         fX0density = 0.425;                    << 312          fX0density = 0.425; fX1density = 2.0; fMdensity = 5.949;
301         fX1density = 2.0;                      << 
302         fMdensity = 5.949;                     << 
303       }                                           313       }
304     }                                          << 314     } else {
305     else {                                     << 
306       //                                          315       //
307       // gases                                    316       // gases
308       //                                          317       //
309       fMdensity = 3.;                             318       fMdensity = 3.;
310       fX1density = 4.0;                           319       fX1density = 4.0;
311                                                << 320       if(fCdensity < 10.) {
312       if (fCdensity <= 10.) {                  << 321   fX0density = 1.6; 
313         fX0density = 1.6;                      << 322       } else if(fCdensity < 11.5) { 
314       }                                        << 323   fX0density = 1.6 + 0.2*(fCdensity - 10.); 
315       else if (fCdensity <= 10.5) {            << 324       } else if(fCdensity < 12.25) { 
316         fX0density = 1.7;                      << 325   fX0density = 1.9 + (fCdensity - 11.5)/7.5; 
317       }                                        << 326       } else if(fCdensity < 13.804) { 
318       else if (fCdensity <= 11.0) {            << 327   fX0density = 2.0; 
319         fX0density = 1.8;                      << 328   fX1density = 4.0 + (fCdensity - 12.25)/1.554;
320       }                                        << 329       } else {
321       else if (fCdensity <= 11.5) {            << 330   fX0density = 0.326*fCdensity-2.5; fX1density = 5.0; 
322         fX0density = 1.9;                      << 
323       }                                           331       }
324       else if (fCdensity <= 12.25) {           << 332       
325         fX0density = 2.0;                      << 333       //special: Hydrogen
326       }                                        << 
327       else if (fCdensity <= 13.804) {          << 
328         fX0density = 2.0;                      << 
329         fX1density = 5.0;                      << 
330       }                                        << 
331       else {                                   << 
332         fX0density = 0.326 * fCdensity - 2.5;  << 
333         fX1density = 5.0;                      << 
334       }                                        << 
335                                                << 
336       // special: Hydrogen                     << 
337       if (1 == nelm && 1 == Z0) {                 334       if (1 == nelm && 1 == Z0) {
338         fX0density = 1.837;                    << 335          fX0density = 1.837; fX1density = 3.0; fMdensity = 4.754;
339         fX1density = 3.0;                      << 
340         fMdensity = 4.754;                     << 
341       }                                           336       }
342                                                << 337       
343       // special: Helium                       << 338       //special: Helium
344       if (1 == nelm && 2 == Z0) {                 339       if (1 == nelm && 2 == Z0) {
345         fX0density = 2.191;                    << 340          fX0density = 2.191; fX1density = 3.0; fMdensity = 3.297;
346         fX1density = 3.0;                      << 
347         fMdensity = 3.297;                     << 
348       }                                           341       }
349     }                                             342     }
350   }                                               343   }
351                                                   344 
352   // change parameters if the gas is not in ST    345   // change parameters if the gas is not in STP.
353   // For the correction the density(STP) is ne << 346   // For the correction the density(STP) is needed. 
354   // Density(STP) is calculated here :         << 347   // Density(STP) is calculated here : 
355                                                << 348       
356   if (State == kStateGas) {                    << 349   if (State == kStateGas) { 
                                                   >> 350     G4double Density  = fMaterial->GetDensity();
357     G4double Pressure = fMaterial->GetPressure    351     G4double Pressure = fMaterial->GetPressure();
358     G4double Temp = fMaterial->GetTemperature( << 352     G4double Temp     = fMaterial->GetTemperature();
359                                                << 353       
360     G4double DensitySTP = density * STP_Pressu << 354     G4double DensitySTP = Density*STP_Pressure*Temp/(Pressure*NTP_Temperature);
361                                                << 355 
362     G4double ParCorr = G4Log(density / Density << 356     G4double ParCorr = G4Log(Density/DensitySTP);
363                                                << 357   
364     fCdensity -= ParCorr;                      << 358     fCdensity  -= ParCorr;
365     fX0density -= ParCorr / twoln10;           << 359     fX0density -= ParCorr/twoln10;
366     fX1density -= ParCorr / twoln10;           << 360     fX1density -= ParCorr/twoln10;
367   }                                            << 361   }
368                                                << 362 
369   // fAdensity parameter can be fixed for not  << 363   // fAdensity parameter can be fixed for not conductive materials 
370   if (0.0 == fD0density) {                     << 364   if(0.0 == fD0density) {
371     G4double Xa = fCdensity / twoln10;         << 365     G4double Xa = fCdensity/twoln10;
372     fAdensity = twoln10 * (Xa - fX0density) /  << 366     fAdensity = twoln10*(Xa-fX0density)
373   }                                            << 367       /std::pow((fX1density-fX0density),fMdensity);
                                                   >> 368   }
                                                   >> 369   /*  
                                                   >> 370   G4cout << "G4IonisParamMat: density effect data for <" 
                                                   >> 371          << fMaterial->GetName() 
                                                   >> 372    << "> " << G4endl;
                                                   >> 373   G4cout << "Eplasma(eV)= " << fPlasmaEnergy/eV
                                                   >> 374    << " rho= " << fAdjustmentFactor
                                                   >> 375    << " -C= " << fCdensity 
                                                   >> 376    << " x0= " << fX0density
                                                   >> 377    << " x1= " << fX1density
                                                   >> 378    << " a= " << fAdensity
                                                   >> 379    << " m= " << fMdensity
                                                   >> 380    << G4endl;
                                                   >> 381   */
374 }                                                 382 }
375                                                   383 
376 //....oooOO0OOooo........oooOO0OOooo........oo    384 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
377                                                   385 
378 void G4IonisParamMat::ComputeFluctModel()         386 void G4IonisParamMat::ComputeFluctModel()
379 {                                                 387 {
380   // compute parameters for the energy loss fl    388   // compute parameters for the energy loss fluctuation model
381   // needs an 'effective Z'                    << 389   // needs an 'effective Z' 
382   G4double Zeff = 0.;                             390   G4double Zeff = 0.;
383   for (std::size_t i = 0; i < fMaterial->GetNu << 391   for (size_t i=0; i<fMaterial->GetNumberOfElements(); ++i) {
384     Zeff += (fMaterial->GetFractionVector())[i << 392      Zeff += (fMaterial->GetFractionVector())[i]
385   }                                            << 393              *((*(fMaterial->GetElementVector()))[i]->GetZ());
386   if (Zeff > 2.1) {                            << 394   }
387     fF2fluct = 2.0 / Zeff;                     << 395   fF2fluct = (Zeff > 2.) ? 2./Zeff : 0.0; 
388     fF1fluct = 1. - fF2fluct;                  << 396 
389     fEnergy2fluct = 10. * Zeff * Zeff * CLHEP: << 397   fF1fluct         = 1. - fF2fluct;
390     fLogEnergy2fluct = G4Log(fEnergy2fluct);   << 398   fEnergy2fluct    = 10.*Zeff*Zeff*eV;
391     fLogEnergy1fluct = (fLogMeanExcEnergy - fF << 399   fLogEnergy2fluct = G4Log(fEnergy2fluct);
392   } else if (Zeff > 1.1) {                     << 400   fLogEnergy1fluct = (fLogMeanExcEnergy - fF2fluct*fLogEnergy2fluct)
393     fF2fluct = 0.0;                            << 401                      /fF1fluct;
394     fF1fluct = 1.0;                            << 402   fEnergy1fluct    = G4Exp(fLogEnergy1fluct);
395     fEnergy2fluct = 40. * CLHEP::eV;           << 403   fEnergy0fluct    = 10.*eV;
396     fLogEnergy2fluct = G4Log(fEnergy2fluct);   << 
397     fLogEnergy1fluct = fLogMeanExcEnergy;      << 
398   } else {                                     << 
399     fF2fluct = 0.0;                            << 
400     fF1fluct = 1.0;                            << 
401     fEnergy2fluct = 10. * CLHEP::eV;           << 
402     fLogEnergy2fluct = G4Log(fEnergy2fluct);   << 
403     fLogEnergy1fluct = fLogMeanExcEnergy;      << 
404   }                                            << 
405   fEnergy1fluct = G4Exp(fLogEnergy1fluct);     << 
406   fEnergy0fluct = 10. * CLHEP::eV;             << 
407   fRateionexcfluct = 0.4;                         404   fRateionexcfluct = 0.4;
408 }                                                 405 }
409                                                   406 
410 //....oooOO0OOooo........oooOO0OOooo........oo    407 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
411                                                   408 
412 void G4IonisParamMat::ComputeIonParameters()      409 void G4IonisParamMat::ComputeIonParameters()
413 {                                                 410 {
414   // get elements in the actual material,         411   // get elements in the actual material,
415   const G4ElementVector* theElementVector = fM << 412   const G4ElementVector* theElementVector = fMaterial->GetElementVector() ;
416   const G4double* theAtomicNumDensityVector =  << 413   const G4double* theAtomicNumDensityVector =
417   const auto NumberOfElements = (G4int)fMateri << 414                          fMaterial->GetAtomicNumDensityVector() ;
                                                   >> 415   const G4int NumberOfElements = fMaterial->GetNumberOfElements() ;
418                                                   416 
419   //  loop for the elements in the material       417   //  loop for the elements in the material
420   //  to find out average values Z, vF, lF        418   //  to find out average values Z, vF, lF
421   G4double z(0.0), vF(0.0), lF(0.0), a23(0.0);    419   G4double z(0.0), vF(0.0), lF(0.0), a23(0.0);
422                                                   420 
423   G4Pow* g4pow = G4Pow::GetInstance();            421   G4Pow* g4pow = G4Pow::GetInstance();
424   if (1 == NumberOfElements) {                 << 422   if( 1 == NumberOfElements ) {
425     const G4Element* element = (*theElementVec    423     const G4Element* element = (*theElementVector)[0];
426     z = element->GetZ();                          424     z = element->GetZ();
427     vF = element->GetIonisation()->GetFermiVel << 425     vF= element->GetIonisation()->GetFermiVelocity();
428     lF = element->GetIonisation()->GetLFactor( << 426     lF= element->GetIonisation()->GetLFactor();
429     a23 = 1.0 / g4pow->A23(element->GetN());   << 427     a23 = 1.0/g4pow->A23(element->GetN());
430   }                                            << 428 
431   else {                                       << 429   } else {
432     G4double norm(0.0);                           430     G4double norm(0.0);
433     for (G4int iel = 0; iel < NumberOfElements << 431     for (G4int iel=0; iel<NumberOfElements; ++iel) {
434       const G4Element* element = (*theElementV    432       const G4Element* element = (*theElementVector)[iel];
435       const G4double weight = theAtomicNumDens    433       const G4double weight = theAtomicNumDensityVector[iel];
436       norm += weight;                             434       norm += weight;
437       z += element->GetZ() * weight;           << 435       z    += element->GetZ() * weight;
438       vF += element->GetIonisation()->GetFermi << 436       vF   += element->GetIonisation()->GetFermiVelocity() * weight;
439       lF += element->GetIonisation()->GetLFact << 437       lF   += element->GetIonisation()->GetLFactor() * weight;
440       a23 += weight / g4pow->A23(element->GetN << 438       a23  += weight/g4pow->A23(element->GetN());
441     }                                          << 439     }
442     if (norm > 0.0) { norm = 1.0/norm; }       << 440     z   /= norm;
443     z *= norm;                                 << 441     vF  /= norm;
444     vF *= norm;                                << 442     lF  /= norm;
445     lF *= norm;                                << 443     a23 /= norm;
446     a23 *= norm;                               << 444   }  
447   }                                            << 445   fZeff        = z;
448   fZeff = z;                                   << 446   fLfactor     = lF;
449   fLfactor = lF;                               << 447   fFermiEnergy = 25.*keV*vF*vF;
450   fFermiEnergy = 25. * CLHEP::keV * vF * vF;   << 448   fInvA23      = a23;
451   fInvA23 = a23;                               << 
452 }                                                 449 }
453                                                   450 
454 //....oooOO0OOooo........oooOO0OOooo........oo    451 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
455                                                   452 
456 void G4IonisParamMat::SetMeanExcitationEnergy(    453 void G4IonisParamMat::SetMeanExcitationEnergy(G4double value)
457 {                                                 454 {
458   if (value == fMeanExcitationEnergy || value  << 455   if(value == fMeanExcitationEnergy || value <= 0.0) { return; }
459     return;                                    << 
460   }                                            << 
461   if (G4NistManager::Instance()->GetVerbose()     456   if (G4NistManager::Instance()->GetVerbose() > 1) {
462     G4cout << "G4Material: Mean excitation ene << 457     G4cout << "G4Material: Mean excitation energy is changed for "
463            << " Iold= " << fMeanExcitationEner << 458            << fMaterial->GetName()
                                                   >> 459            << " Iold= " << fMeanExcitationEnergy/eV
                                                   >> 460            << "eV; Inew= " << value/eV << " eV;"
464            << G4endl;                             461            << G4endl;
465   }                                               462   }
466                                                << 463   
467   fMeanExcitationEnergy = value;                  464   fMeanExcitationEnergy = value;
468                                                   465 
469   // add corrections to density effect            466   // add corrections to density effect
470   G4double newlog = G4Log(value);                 467   G4double newlog = G4Log(value);
471   G4double corr = 2 * (newlog - fLogMeanExcEne << 468   G4double corr = 2*(newlog - fLogMeanExcEnergy);
472   fCdensity += corr;                           << 469   fCdensity  += corr;
473   fX0density += corr / twoln10;                << 470   fX0density += corr/twoln10;
474   fX1density += corr / twoln10;                << 471   fX1density += corr/twoln10;
475                                                   472 
476   // recompute parameters of fluctuation model    473   // recompute parameters of fluctuation model
477   fLogMeanExcEnergy = newlog;                     474   fLogMeanExcEnergy = newlog;
478   ComputeFluctModel();                            475   ComputeFluctModel();
479 }                                                 476 }
480                                                   477 
481 //....oooOO0OOooo........oooOO0OOooo........oo    478 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
482                                                   479 
483 void G4IonisParamMat::SetDensityEffectParamete << 480 void G4IonisParamMat::SetDensityEffectParameters(G4double cd, G4double md, 
484   G4double cd, G4double md, G4double ad, G4dou << 481              G4double ad, G4double x0, 
485 {                                              << 482              G4double x1, G4double d0)
486   // no check on consistence of user parameter << 483 {
487   G4AutoLock l(&ionisMutex);                   << 484   // no check on consistence of user parameters  
488   fCdensity = cd;                              << 485 #ifdef G4MULTITHREADED
489   fMdensity = md;                              << 486   G4MUTEXLOCK(&ionisMutex);
490   fAdensity = ad;                              << 487 #endif
                                                   >> 488   fCdensity  = cd; 
                                                   >> 489   fMdensity  = md;
                                                   >> 490   fAdensity  = ad;
491   fX0density = x0;                                491   fX0density = x0;
492   fX1density = x1;                                492   fX1density = x1;
493   fD0density = d0;                                493   fD0density = d0;
494   l.unlock();                                  << 494 #ifdef G4MULTITHREADED
                                                   >> 495   G4MUTEXUNLOCK(&ionisMutex);
                                                   >> 496 #endif
495 }                                                 497 }
496                                                   498 
497 //....oooOO0OOooo........oooOO0OOooo........oo    499 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
498                                                   500 
499 void G4IonisParamMat::SetDensityEffectParamete    501 void G4IonisParamMat::SetDensityEffectParameters(const G4Material* bmat)
500 {                                                 502 {
501   G4AutoLock l(&ionisMutex);                   << 503 #ifdef G4MULTITHREADED
                                                   >> 504   G4MUTEXLOCK(&ionisMutex);
                                                   >> 505 #endif
502   const G4IonisParamMat* ipm = bmat->GetIonisa    506   const G4IonisParamMat* ipm = bmat->GetIonisation();
503   fCdensity = ipm->GetCdensity();              << 507   fCdensity  = ipm->GetCdensity(); 
504   fMdensity = ipm->GetMdensity();              << 508   fMdensity  = ipm->GetMdensity();
505   fAdensity = ipm->GetAdensity();              << 509   fAdensity  = ipm->GetAdensity();
506   fX0density = ipm->GetX0density();               510   fX0density = ipm->GetX0density();
507   fX1density = ipm->GetX1density();               511   fX1density = ipm->GetX1density();
508   fD0density = ipm->GetD0density();               512   fD0density = ipm->GetD0density();
509                                                   513 
510   G4double corr = G4Log(bmat->GetDensity() / f << 514   G4double corr = G4Log(bmat->GetDensity()/fMaterial->GetDensity());
511   fCdensity += corr;                           << 515   fCdensity  += corr;
512   fX0density += corr / twoln10;                << 516   fX0density += corr/twoln10;
513   fX1density += corr / twoln10;                << 517   fX1density += corr/twoln10;
514   l.unlock();                                  << 518 #ifdef G4MULTITHREADED
                                                   >> 519   G4MUTEXUNLOCK(&ionisMutex);
                                                   >> 520 #endif
515 }                                                 521 }
516                                                   522 
517 //....oooOO0OOooo........oooOO0OOooo........oo    523 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
518                                                   524 
519 void G4IonisParamMat::ComputeDensityEffectOnFl    525 void G4IonisParamMat::ComputeDensityEffectOnFly(G4bool val)
520 {                                                 526 {
521   if (val) {                                   << 527   if(val) {
522     if (nullptr == fDensityEffectCalc) {       << 528     if(!fDensityEffectCalc) { 
523       G4int n = 0;                                529       G4int n = 0;
524       for (std::size_t i = 0; i < fMaterial->G << 530       for(size_t i=0; i<fMaterial->GetNumberOfElements(); ++i) {
525         const G4int Z = fMaterial->GetElement( << 531   const G4int Z = fMaterial->GetElement(i)->GetZasInt();
526         n += G4AtomicShells::GetNumberOfShells    532         n += G4AtomicShells::GetNumberOfShells(Z);
527       }                                           533       }
528       // The last level is the conduction leve    534       // The last level is the conduction level.  If this is *not* a conductor,
529       // make a dummy conductor level with zer    535       // make a dummy conductor level with zero electrons in it.
530       fDensityEffectCalc = new G4DensityEffect    536       fDensityEffectCalc = new G4DensityEffectCalculator(fMaterial, n);
531     }                                             537     }
532   }                                            << 538   } else {
533   else {                                       << 
534     delete fDensityEffectCalc;                    539     delete fDensityEffectCalc;
535     fDensityEffectCalc = nullptr;                 540     fDensityEffectCalc = nullptr;
536   }                                               541   }
537 }                                                 542 }
538                                                   543 
539 //....oooOO0OOooo........oooOO0OOooo........oo    544 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
540                                                   545 
541 G4double G4IonisParamMat::FindMeanExcitationEn    546 G4double G4IonisParamMat::FindMeanExcitationEnergy(const G4Material* mat) const
542 {                                                 547 {
543   G4double res = 0.0;                             548   G4double res = 0.0;
544   // data from density effect data                549   // data from density effect data
545   if (fDensityData != nullptr) {               << 550   if(fDensityData) {
546     G4int idx = fDensityData->GetIndex(mat->Ge    551     G4int idx = fDensityData->GetIndex(mat->GetName());
547     if (idx >= 0) {                            << 552     if(idx >= 0) {
548       res = fDensityData->GetMeanIonisationPot    553       res = fDensityData->GetMeanIonisationPotential(idx);
549     }                                             554     }
550   }                                               555   }
551                                                   556 
552   // The data on mean excitation energy for co    557   // The data on mean excitation energy for compaunds
553   // from "Stopping Powers for Electrons and P    558   // from "Stopping Powers for Electrons and Positrons"
554   // ICRU Report N#37, 1984  (energy in eV)       559   // ICRU Report N#37, 1984  (energy in eV)
555   // this value overwrites Density effect data    560   // this value overwrites Density effect data
556   G4String chFormula = mat->GetChemicalFormula    561   G4String chFormula = mat->GetChemicalFormula();
557   if (! chFormula.empty()) {                   << 562   if(chFormula != "") { 
558     static const size_t numberOfMolecula = 54; << 563 
559     // clang-format off                        << 564     static const size_t numberOfMolecula = 54; 
560     static const G4String name[numberOfMolecul    565     static const G4String name[numberOfMolecula] = {
561       // gas 0 - 12                               566       // gas 0 - 12
562       "NH_3",       "C_4H_10",     "CO_2",        567       "NH_3",       "C_4H_10",     "CO_2",      "C_2H_6",      "C_7H_16-Gas",
563       // "G4_AMMONIA", "G4_BUTANE","G4_CARBON_    568       // "G4_AMMONIA", "G4_BUTANE","G4_CARBON_DIOXIDE","G4_ETHANE", "G4_N-HEPTANE"
564       "C_6H_14-Gas",   "CH_4",     "NO",          569       "C_6H_14-Gas",   "CH_4",     "NO",        "N_2O",        "C_8H_18-Gas",      
565       // "G4_N-HEXANE" , "G4_METHANE", "x", "G    570       // "G4_N-HEXANE" , "G4_METHANE", "x", "G4_NITROUS_OXIDE", "G4_OCTANE"
566       "C_5H_12-Gas",   "C_3H_8",   "H_2O-Gas",    571       "C_5H_12-Gas",   "C_3H_8",   "H_2O-Gas",                        
567       // "G4_N-PENTANE", "G4_PROPANE", "G4_WAT    572       // "G4_N-PENTANE", "G4_PROPANE", "G4_WATER_VAPOR"
568                                                   573 
569       // liquid 13 - 39                           574       // liquid 13 - 39
570       "C_3H_6O",    "C_6H_5NH_2",  "C_6H_6",      575       "C_3H_6O",    "C_6H_5NH_2",  "C_6H_6",    "C_4H_9OH",    "CCl_4", 
571       //"G4_ACETONE","G4_ANILINE","G4_BENZENE" << 576       //"G4_ACETONE","G4_ANILINE","G4_BENZENE","G4_N-BUTYL_ALCOHOL","G4_CARBON_TETRACHLORIDE"  
572       "C_6H_5Cl",   "CHCl_3",      "C_6H_12",     577       "C_6H_5Cl",   "CHCl_3",      "C_6H_12",   "C_6H_4Cl_2",  "C_4Cl_2H_8O",
573       //"G4_CHLOROBENZENE","G4_CHLOROFORM","G4    578       //"G4_CHLOROBENZENE","G4_CHLOROFORM","G4_CYCLOHEXANE","G4_1,2-DICHLOROBENZENE",
574       //"G4_DICHLORODIETHYL_ETHER"                579       //"G4_DICHLORODIETHYL_ETHER"
575       "C_2Cl_2H_4", "(C_2H_5)_2O", "C_2H_5OH",    580       "C_2Cl_2H_4", "(C_2H_5)_2O", "C_2H_5OH",  "C_3H_5(OH)_3","C_7H_16",
576       //"G4_1,2-DICHLOROETHANE","G4_DIETHYL_ET    581       //"G4_1,2-DICHLOROETHANE","G4_DIETHYL_ETHER","G4_ETHYL_ALCOHOL","G4_GLYCEROL","G4_N-HEPTANE"
577       "C_6H_14",    "CH_3OH",      "C_6H_5NO_2    582       "C_6H_14",    "CH_3OH",      "C_6H_5NO_2","C_5H_12",     "C_3H_7OH",
578       //"G4_N-HEXANE","G4_METHANOL","G4_NITROB << 583       //"G4_N-HEXANE","G4_METHANOL","G4_NITROBENZENE","G4_N-PENTANE","G4_N-PROPYL_ALCOHOL",  
579       "C_5H_5N",    "C_8H_8",      "C_2Cl_4",     584       "C_5H_5N",    "C_8H_8",      "C_2Cl_4",   "C_7H_8",      "C_2Cl_3H",
580       //"G4_PYRIDINE","G4_POLYSTYRENE","G4_TET    585       //"G4_PYRIDINE","G4_POLYSTYRENE","G4_TETRACHLOROETHYLENE","G4_TOLUENE","G4_TRICHLOROETHYLENE"
581       "H_2O",       "C_8H_10",                    586       "H_2O",       "C_8H_10",             
582       // "G4_WATER", "G4_XYLENE"                  587       // "G4_WATER", "G4_XYLENE"
583                                                   588 
584       // solid 40 - 53                            589       // solid 40 - 53
585       "C_5H_5N_5",  "C_5H_5N_5O",  "(C_6H_11NO    590       "C_5H_5N_5",  "C_5H_5N_5O",  "(C_6H_11NO)-nylon",  "C_25H_52",
586       // "G4_ADENINE", "G4_GUANINE", "G4_NYLON    591       // "G4_ADENINE", "G4_GUANINE", "G4_NYLON-6-6", "G4_PARAFFIN"
587       "(C_2H_4)-Polyethylene",     "(C_5H_8O_2    592       "(C_2H_4)-Polyethylene",     "(C_5H_8O_2)-Polymethil_Methacrylate",
588       // "G4_ETHYLENE", "G4_PLEXIGLASS"           593       // "G4_ETHYLENE", "G4_PLEXIGLASS"
589       "(C_8H_8)-Polystyrene",      "A-150-tiss    594       "(C_8H_8)-Polystyrene",      "A-150-tissue",       "Al_2O_3",  "CaF_2",
590       // "G4_POLYSTYRENE", "G4_A-150_TISSUE",     595       // "G4_POLYSTYRENE", "G4_A-150_TISSUE", "G4_ALUMINUM_OXIDE", "G4_CALCIUM_FLUORIDE"
591       "LiF",        "Photo_Emulsion",  "(C_2F_    596       "LiF",        "Photo_Emulsion",  "(C_2F_4)-Teflon",  "SiO_2"
592       // "G4_LITHIUM_FLUORIDE", "G4_PHOTO_EMUL    597       // "G4_LITHIUM_FLUORIDE", "G4_PHOTO_EMULSION", "G4_TEFLON", "G4_SILICON_DIOXIDE"
593     } ;                                           598     } ;
594                                                << 599     
595     static const G4double meanExcitation[numbe    600     static const G4double meanExcitation[numberOfMolecula] = {
596                                                   601 
597       53.7,   48.3,  85.0,  45.4,  49.2,          602       53.7,   48.3,  85.0,  45.4,  49.2,
598       49.1,   41.7,  87.8,  84.9,  49.5,          603       49.1,   41.7,  87.8,  84.9,  49.5,
599       48.2,   47.1,  71.6,                        604       48.2,   47.1,  71.6,
600                                                   605 
601       64.2,   66.2,  63.4,  59.9,  166.3,         606       64.2,   66.2,  63.4,  59.9,  166.3,
602       89.1,  156.0,  56.4, 106.5,  103.3,         607       89.1,  156.0,  56.4, 106.5,  103.3, 
603       111.9,   60.0,  62.9,  72.6,   54.4,        608       111.9,   60.0,  62.9,  72.6,   54.4,  
604       54.0,  67.6,   75.8,  53.6,   61.1,         609       54.0,  67.6,   75.8,  53.6,   61.1,  
605       66.2,  64.0,  159.2,  62.5,  148.1,         610       66.2,  64.0,  159.2,  62.5,  148.1,  
606       75.0,  61.8,                                611       75.0,  61.8,
607                                                   612 
608       71.4,  75.0,   63.9,  48.3,   57.4,         613       71.4,  75.0,   63.9,  48.3,   57.4,
609       74.0,  68.7,   65.1, 145.2,  166.,          614       74.0,  68.7,   65.1, 145.2,  166.,
610       94.0, 331.0,   99.1, 139.2                  615       94.0, 331.0,   99.1, 139.2 
611     };                                            616     };
612     // clang-format on                         << 
613                                                   617 
614     for (std::size_t i = 0; i < numberOfMolecu << 618     for(size_t i=0; i<numberOfMolecula; i++) {
615       if (chFormula == name[i]) {              << 619       if(chFormula == name[i]) {
616         res = meanExcitation[i] * CLHEP::eV;   << 620   res = meanExcitation[i]*eV;
617         break;                                 << 621   break;
618       }                                           622       }
619     }                                             623     }
620   }                                               624   }
621   return res;                                     625   return res;
622 }                                                 626 }
                                                   >> 627 
                                                   >> 628 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
                                                   >> 629 
623                                                   630