In Physics Reference Manual (p.159) and original source of R.M. Sternheimer and R.F. Peierls (General expression for the density effect for the ionization lossof charged particles. Physical Review B, 3(11):3681–3692, jun 1971. URL: https://doi.org/10.1103/PhysRevB.3.3681, doi:10.1103/physrevb.3.3681) parameters x0 and x1 are discrete: for C < 10 x0 = 1.6, x1 = 4 for C in [10.0, 10.5[ x0 = 1.7, x1 = 4 for C in [10.5, 11.0[ x0 = 1.8, x1 = 4 for C in [11.0, 11.5[ x0 = 1.9, x1 = 4 ... so on ... But in the function G4IonisParamMat::ComputeDensityEffectParameters (lines 322-328 of file G4IonisParamMat.cc) linear interpolation for these parameters is used thus making the parameters wrong (inconsistent with original source) inside given intervals. This can be illustrated by the following code. Its output shows that (epsilon << 1): for C = 10.0 + epsilon parameter x0 ~ 1.6 (instead of 1.7) for C = 10.5 + epsilon parameter x0 ~ 1.7 (instead of 1.8) for C = 11.0 + epsilon parameter x0 ~ 1.8 (instead of 1.9) ... so on ... ********** CODE ********** G4NistManager* man = G4NistManager::Instance(); G4Element* elH = man->FindOrBuildElement("H"); G4Element* elO = man->FindOrBuildElement("O"); constexpr G4double gcc = CLHEP::g / CLHEP::cm3; G4double rho0 = 0.1 * gcc; // Cdensity = const - log(rho) std::stringstream name0; name0 << "H2O" << "(" << rho0 / gcc << ")"; G4Material* mat0 = new G4Material(name0.str(), rho0, 2, kStateGas); mat0->AddElement(elH, 2); mat0->AddElement(elO, 1); G4double C0 = mat0->GetIonisation()->GetCdensity(); constexpr G4double eps = 1.0e-3; std::vector<G4double> Cs{10.0, 10.5, 11.0, 11.5, 12.25, 13.804}; for(G4double C : Cs) { G4double C1 = C - eps; G4double C2 = C + eps; G4double rho1 = rho0 * std::exp(C0 - C1); G4double rho2 = rho0 * std::exp(C0 - C2); std::stringstream name1; name1 << "H2O" << "(" << rho1 / gcc << ")"; G4Material* mat1 = new G4Material(name1.str(), rho1, 2, kStateGas); mat1->AddElement(elH, 2); mat1->AddElement(elO, 1); std::stringstream name2; name2 << "H2O" << "(" << rho2 / gcc << ")"; G4Material* mat2 = new G4Material(name2.str(), rho2, 2, kStateGas); mat2->AddElement(elH, 2); mat2->AddElement(elO, 1); G4cout << C << "+"; G4cout << "(" << mat1->GetIonisation()->GetCdensity() - C << ")" << "\t"; G4cout << mat1->GetIonisation()->GetX0density() << G4endl; G4cout << C << "+"; G4cout << "(" << mat2->GetIonisation()->GetCdensity() - C << ")" << "\t"; G4cout << mat2->GetIonisation()->GetX0density() << G4endl; G4cout << G4endl; } ********** OUTPUT ********** 10+(-0.001) 1.6 10+(0.001) 1.6002 10.5+(-0.001) 1.6998 10.5+(0.001) 1.7002 11+(-0.001) 1.7998 11+(0.001) 1.8002 11.5+(-0.001) 1.8998 11.5+(0.001) 1.90013 12.25+(-0.001) 1.99987 12.25+(0.001) 2 13.804+(-0.001) 2 13.804+(0.001) 2.00043
Hello, we decided return back version of this place from Geant4 9.5,following by your proposal. The effect is small but we agree that better to be fully consistent with the original paper. The fix will be available with the next public release and will be backported to the next patch to 10.7. VI