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 3.1)


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