Problem 2334

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