//

// ********************************************************************

// * License and Disclaimer                                           *

// *                                                                  *

// * The  Geant4 software  is  copyright of the Copyright Holders  of *

// * the Geant4 Collaboration.  It is provided  under  the terms  and *

// * conditions of the Geant4 Software License,  included in the file *

// * LICENSE and available at  http://cern.ch/geant4/license .  These *

// * include a list of copyright holders.                             *

// *                                                                  *

// * Neither the authors of this software system, nor their employing *

// * institutes,nor the agencies providing financial support for this *

// * work  make  any representation or  warranty, express or implied, *

// * regarding  this  software system or assume any liability for its *

// * use.  Please see the license in the file  LICENSE  and URL above *

// * for the full disclaimer and the limitation of liability.         *

// *                                                                  *

// * This  code  implementation is the result of  the  scientific and *

// * technical work of the GEANT4 collaboration.                      *

// * By using,  copying,  modifying or  distributing the software (or *

// * any work based  on the software)  you  agree  to acknowledge its *

// * use  in  resulting  scientific  publications,  and indicate your *

// * acceptance of all terms of the Geant4 Software license.          *

// ********************************************************************

//

/*

 * Interface to calculation of the Fermi density effect as per the method

 * described in:

 *

 *   R. M. Sternheimer, M. J. Berger, and S. M. Seltzer. Density

 *   effect for the ionization loss of charged particles in various sub-

 *   stances. Atom. Data Nucl. Data Tabl., 30:261, 1984.

 *

 * Which (among other Sternheimer references) builds on:

 *

 *   R. M. Sternheimer. The density effect for ionization loss in

 *   materials. Phys. Rev., 88:851­859, 1952.

 *

 * The returned values of delta are directly from the Sternheimer calculation,

 * and not Sternheimer's popular three-part approximate parameterization

 * introduced in the same paper.

 *

 * Author: Matthew Strait <straitm@umn.edu> 2019

 */

#ifndef G4DensityEffectCalc_HH

#define G4DensityEffectCalc_HH

/*

 * Holds the user input and intermediate steps for calculation of

 * the Fermi density effect parameter delta, a term in the Bethe

 * energy loss formula.

 */

struct G4DensityEffectCalcData{

  /***** User set *****/

  // Number of energy levels.  If a single element, this is the number

  // of subshells.  If several elements, this is the sum of the number

  // of subshells.  In principle, could include levels for molecular

  // orbitals or other non-atomic states.  The last level is always

  // the conduction band.  If the material is an insulator, set the

  // oscillator strength for that level to zero and the energy to

  // any value.

  int nlev;

  // Sternheimer's "oscillator strengths", which are simply the fraction

  // of electrons in a given energy level.  For a single element, this is

  // the fraction of electrons in a subshell.  For a compound or mixture,

  // it is weighted by the number fraction of electrons contributed by

  // each element, e.g. for water, oxygen's electrons are given 8/10 of the

  // weight.

  double * sternf;

  // Energy levels.  Can be found for free atoms in, e.g., T. A. Carlson.

  // Photoelectron and Auger Spectroscopy. Plenum Press, New York and London,

  // 1985. Available in a convenient form in G4AtomicShells.cc.

  //

  // Sternheimer 1984 implies that the energy level for conduction electrons

  // (the final element of this array) should be set to zero, although the

  // computation could be run with other values.

  double * levE;

  // The plasma energy of the material in eV, which is simply

  // 28.816 sqrt(density Z/A), with density in g/cc.

  double plasmaE;

  // The mean excitation energy of the material in eV, i.e. the 'I' in the

  // Bethe energy loss formula.

  double meanexcite;

  /***** Results of intermediate calculations *****/

  // The Sternheimer parameters l_i which appear in Sternheimer 1984 eq(1).

  double * sternl;

  // The adjusted energy levels, as found using Sternheimer 1984 eq(8).

  double * sternEbar;

  /***** Packaged for convenience *****/

  // The Sternheimer 'x' defined as log10(p/m) == log10(beta*gamma).

  // This is packaged with the rest of the data for convenience.

  // Any value set here is overwritten when the user calls DoFermiDeltaCalc().

  double sternx;

};

/**

 * Given a material defined in 'par' with a plasma energy, mean excitation

 * energy, and set of atomic energy levels ("oscillator frequencies") with

 * occupation fractions ("oscillation strengths"), solve for the Sternheimer

 * adjustment factor (Sternheimer 1984 eq 8) and record (into 'par') the values

 * of the adjusted oscillator frequencies and Sternheimer constants l_i.

 * After doing this, 'par' is ready for a calculation of delta for an

 * arbitrary particle energy.  Returns true on success, false on failure.

 */

bool SetupFermiDeltaCalc(G4DensityEffectCalcData * par);

/**

 * Given that SetupFermiDeltaCalc() has already been called, calculate

 * the Fermi density effect delta for a given Sternheimer x, defined as

 * log10(beta * gamma).  This can be called any number of times without

 * rerunning SetupFermiDeltaCalc().

 *

 * Physically, delta is always non-negative.  Returns -1 on failure.

 */

double DoFermiDeltaCalc(G4DensityEffectCalcData * par,

                        const double sternx);

#endif

#include "G4DensityEffectCalc.hh"

  // Stores atomic data necessary for computing the Sternheimer exact

  // density effect, as well as intermediate stages of the computation.

  G4DensityEffectCalcData * fCalcDensity;

  // If the computation of the Sternheimer exact density effect fails,

  // this flag is set to prevent wasting time on repeated failures.

  G4bool fCalcDensityFailure = false;

			G4DensityEffectCalc.hh

			G4DensityEffectCalc.cc

//

// ********************************************************************

// * License and Disclaimer                                           *

// *                                                                  *

// * The  Geant4 software  is  copyright of the Copyright Holders  of *

// * the Geant4 Collaboration.  It is provided  under  the terms  and *

// * conditions of the Geant4 Software License,  included in the file *

// * LICENSE and available at  http://cern.ch/geant4/license .  These *

// * include a list of copyright holders.                             *

// *                                                                  *

// * Neither the authors of this software system, nor their employing *

// * institutes,nor the agencies providing financial support for this *

// * work  make  any representation or  warranty, express or implied, *

// * regarding  this  software system or assume any liability for its *

// * use.  Please see the license in the file  LICENSE  and URL above *

// * for the full disclaimer and the limitation of liability.         *

// *                                                                  *

// * This  code  implementation is the result of  the  scientific and *

// * technical work of the GEANT4 collaboration.                      *

// * By using,  copying,  modifying or  distributing the software (or *

// * any work based  on the software)  you  agree  to acknowledge its *

// * use  in  resulting  scientific  publications,  and indicate your *

// * acceptance of all terms of the Geant4 Software license.          *

// ********************************************************************

//

/*

 * Implements calculation of the Fermi density effect as per the method

 * described in:

 *

 *   R. M. Sternheimer, M. J. Berger, and S. M. Seltzer. Density

 *   effect for the ionization loss of charged particles in various sub-

 *   stances. Atom. Data Nucl. Data Tabl., 30:261, 1984.

 *

 * Which (among other Sternheimer references) builds on:

 *

 *   R. M. Sternheimer. The density effect for ionization loss in

 *   materials. Phys. Rev., 88:851­859, 1952.

 *

 * The returned values of delta are directly from the Sternheimer calculation,

 * and not Sternheimer's popular three-part approximate parameterization

 * introduced in the same paper.

 *

 * Author: Matthew Strait <straitm@umn.edu> 2019

 */

#include "G4ios.hh"

#include "G4DensityEffectCalc.hh"

#include <stdio.h>

#include <stdlib.h>

#include <math.h>

/* Newton's method for finding roots.  Adapted from G4PolynominalSolver, but

 * without the assumption that the input is a polynomial.  Also, here we

 * always expect the roots to be positive, so return -1 as an error value. */

static double newton(const double start,

                     double(*Function)(double, G4DensityEffectCalcData *),

                     double(*Derivative)(double, G4DensityEffectCalcData *),

                     G4DensityEffectCalcData * par)

{

  const int maxIter = 100;

  int nbad = 0, ngood = 0;

  double Lambda = start;

  while(true){

    const double Value = Function(Lambda, par);

    const double Gradient = Derivative(Lambda, par);

    Lambda -= Value/Gradient;

    if(std::fabs((Value/Gradient)/Lambda) <= 1e-12) {

      ngood++;

      if(ngood == 2) return Lambda;

    } else {

      nbad++;

      if(nbad > maxIter || isnan(Value) || isinf(Value)) return -1;

    }

  }

}

/* Return the derivative of the equation used

 * to solve for the Sternheimer parameter rho. */

static double dfrho(double rho, G4DensityEffectCalcData * p)

{

  double ans = 0;

  for(int i = 0; i < p->nlev-1; i++)

    if(p->sternf[i] > 0)

      ans += p->sternf[i] * pow(p->levE[i], 2) * rho /

        (pow(p->levE[i] * rho, 2) + 2./3. * p->sternf[i] * pow(p->plasmaE, 2));

  return ans;

}

/* Return the functional value for the equation used

 * to solve for the Sternheimer parameter rho. */

static double frho(double rho, G4DensityEffectCalcData * p)

{

  double ans = 0;

  for(int i = 0; i < p->nlev-1; i++)

    if(p->sternf[i] > 0)

      ans += p->sternf[i] * log(pow(p->levE[i]*rho, 2) +

        2./3. * p->sternf[i]*pow(p->plasmaE, 2));

  ans *= 0.5; // pulled out of loop for efficiency

  if(p->sternf[p->nlev-1] > 0)

    ans += p->sternf[p->nlev-1] * log(p->plasmaE * sqrt(p->sternf[p->nlev-1]));

  ans -= log(p->meanexcite);

  return ans;

}

/* Return the derivative for the equation used to

 * solve for the Sternheimer parameter l, called 'L' here. */

static double dell(double L, G4DensityEffectCalcData * p)

{

  double ans = 0;

  for(int i = 0; i < p->nlev; i++)

    if(p->sternf[i] > 0 && (p->sternEbar[i] > 0 || L != 0))

      ans += p->sternf[i]/pow(pow(p->sternEbar[i], 2) + L*L, 2);

  ans *= -2*L; // pulled out of the loop for efficiency

  return ans;

}

/* Return the functional value for the equation used to

 * solve for the Sternheimer parameter l, called 'L' here. */

static double ell(double L, G4DensityEffectCalcData * p)

{

  double ans = 0;

  for(int i = 0; i < p->nlev; i++)

    if(p->sternf[i] > 0 && (p->sternEbar[i] > 0 || L != 0))

      ans += p->sternf[i]/(pow(p->sternEbar[i], 2) + L*L);

  ans -= pow(10, - 2 * p->sternx);

  return ans;

}

/**

 * Given the Sternheimer parameter l^2 (called 'sternL' here), and that

 * the l_i and adjusted energies have been found with SetupFermiDeltaCalc(),

 * return the value of delta.  Helper function for DoFermiDeltaCalc().

 */

static double delta_once_solved(G4DensityEffectCalcData * p,

                                const double sternL)

{

  double ans = 0;

  for(int i = 0; i < p->nlev; i++)

    if(p->sternf[i] > 0)

      ans += p->sternf[i] *

        log((pow(p->sternl[i], 2) + pow(sternL, 2))/pow(p->sternl[i], 2));

  ans -= pow(sternL, 2)/(1 + pow(10, 2 * p->sternx));

  return ans;

}

/**

 * Calculate the Sternheimer adjusted energy levels and parameters l_i given

 * the Sternheimer parameter rho.  Helper function for SetupFermiDeltaCalc().

 */

static void calc_Ebarl(G4DensityEffectCalcData * par,

                       const double sternrho)

{

  par->sternl    = (double *)malloc(sizeof(double)*par->nlev);

  par->sternEbar = (double *)malloc(sizeof(double)*par->nlev);

  for(int i = 0; i < par->nlev; i++){

    par->sternEbar[i] = par->levE[i] * sternrho/par->plasmaE;

    par->sternl[i] = i < par->nlev-1?

      sqrt(pow(par->sternEbar[i], 2) + 2./3. * par->sternf[i])

      : sqrt(par->sternf[i]);

  }

}

/**

 * If the "oscillator strengths" par->sternf do not add up to 1,

 * normalize them so that they do.

 */

static void normalize_sternf(G4DensityEffectCalcData * par)

{

  double sum = 0;

  for(int i = 0; i < par->nlev; i++) sum += par->sternf[i];

  for(int i = 0; i < par->nlev; i++) par->sternf[i] /= sum;

}

/**

 * At sufficiently high energy, the density effect takes on a simple

 * form.  Also at sufficiently high energy, we run into numerical problems

 * doing the full calculation.  In that case, return this.

 */

static double delta_limiting_case(G4DensityEffectCalcData * par,

                                  const double sternx)

{

  return 2 * (log(par->plasmaE/par->meanexcite) + log(10.)*sternx) - 1;

}

/********************/

/* Public functions */

/********************/

bool SetupFermiDeltaCalc(G4DensityEffectCalcData * par)

{

  normalize_sternf(par);

  const double sternrho = newton(1.5, frho, dfrho, par);

  // Negative values, and values much larger than unity are non-physical.

  // Values between zero and one are also suspect, but not as clearly wrong.

  if(sternrho <= 0 || sternrho > 100){

    G4cerr <<

    "*********************************************************************\n"

    "* Warning: Could not solve for Sternheimer rho. Probably you have a *\n"

    "* mean ionization energy which is incompatible with your            *\n"

    "* distribution of energy levels.  Falling back to parameterization. *\n"

    "*********************************************************************\n";

    return false;

  }

  calc_Ebarl(par, sternrho);

  return true;

}

double DoFermiDeltaCalc(G4DensityEffectCalcData * par,

                        const double sternx)

{

  // Above beta*gamma of 10^10, the exact treatment is within machine

  // precision of the limiting case, for ordinary solids, at least. The

  // convergence goes up as the density goes down, but even in a pretty

  // hard vacuum it converges by 10^20. Also, it's hard to imagine how

  // this energy is relevant (x = 20 -> 10^19 GeV for muons). So this

  // is mostly not here for physical reasons, but rather to avoid ugly

  // discontinuities in the return value.

  if(sternx > 20) return delta_limiting_case(par, sternx);

  par->sternx = sternx;

  // The derivative of the function we are solving for is strictly

  // negative for positive (physical) values, so if the value at

  // zero is less than zero, it has no solution, and there is no

  // density effect in the Sternheimer "exact" treatment (which is

  // still an approximation).

  if(ell(0, par) <= 0) return 0;

  // Increase initial guess until it converges.

  int startLi;

  for(startLi = -10; startLi < 30; startLi++){

    const double sternL = newton(pow(2, startLi), &ell, &dell, par);

    if(sternL != -1)

      return delta_once_solved(par, sternL);

  }

  return -1; // Signal the caller to use the Sternheimer approximation,

             // because we have been unable to solve the exact form.

}

// 16-01-19, add exact computation of the density effect (M. Strait)

#include "G4DensityEffectCalc.hh"

#include "G4AtomicShells.hh"

  fCalcDensity = nullptr;

  fCalcDensity = nullptr;

void G4IonisParamMat::SetSternheimerExactDensityEffect()

{

  std::vector<G4int> Z;

  std::vector<double> numberfracs;

  const bool isconductor =

       fMaterial->GetMaterialPropertiesTable() != nullptr

   && fMaterial->GetMaterialPropertiesTable()->ConstPropertyExists("conductor")

   && fMaterial->GetMaterialPropertiesTable()->GetConstProperty("conductor");

  if(fCalcDensity != nullptr) delete fCalcDensity;

  fCalcDensity = new G4DensityEffectCalcData;

  memset(fCalcDensity, 0, sizeof(G4DensityEffectCalcData));

  fCalcDensity->nlev = 0;

  for(unsigned int i = 0; i < fMaterial->GetNumberOfElements(); i++){

    Z.push_back(fMaterial->GetElement(i)->GetZ());

    fCalcDensity->nlev += G4AtomicShells::GetNumberOfShells(Z[i]);

  }

  // The last level is the conduction level.  If this is *not* a conductor,

  // make a dummy conductor level with zero electrons in it.

  if(!isconductor) fCalcDensity->nlev++;

  for(unsigned int i = 0; i < fMaterial->GetNumberOfElements(); i++)

    numberfracs.push_back(fMaterial->GetVecNbOfAtomsPerVolume()[i]/

                          fMaterial->GetTotNbOfAtomsPerVolume());

  fCalcDensity->sternf = (double *)malloc(sizeof(double) * fCalcDensity->nlev);

  fCalcDensity->levE   = (double *)malloc(sizeof(double) * fCalcDensity->nlev);

  memset(fCalcDensity->levE, 0, sizeof(double)*fCalcDensity->nlev);

  int sh = 0;

  for(unsigned int j = 0; j < fMaterial->GetNumberOfElements(); j++){

    // The last subshell is considered to contain the conduction

    // electrons. Sternheimer 1984 says "the lowest chemical valance of

    // the element" is used to set the number of conduction electrons.

    // I'm not sure if that means the highest subshell or the whole

    // shell, but in any case, he also says that the choice is arbitrary

    // and offers a possible alternative. This is one of the sources of

    // uncertainty in the model.

    const int nshell = G4AtomicShells::GetNumberOfShells(Z[j]);

    for(int i = 0; i < nshell; i++){

      // For conductors, put *all* top shell electrons into the conduction

      // band, regardless of element.

      const int lev = i < nshell-1 || !isconductor? sh: fCalcDensity->nlev-1;

      fCalcDensity->sternf[lev] += numberfracs[j] *

        double(G4AtomicShells::GetNumberOfElectrons(Z[j], i))/Z[j];

      fCalcDensity->levE[lev] = G4AtomicShells::GetBindingEnergy(Z[j], i)/eV;

      sh++;

    }

  }

  if(!isconductor){

    fCalcDensity->sternf[fCalcDensity->nlev-1] = 0; // No electrons

    fCalcDensity->levE[fCalcDensity->nlev-1] = 1; // dummy value

  }

  else{

    fCalcDensity->levE[fCalcDensity->nlev-1] = 0;

  }

  double sumZ = 0, sumA = 0;

  for(unsigned int i = 0; i < fMaterial->GetNumberOfElements(); i++){

    sumA += numberfracs[i] * fMaterial->GetElement(i)->GetAtomicMassAmu();

    sumZ += numberfracs[i] * fMaterial->GetElement(i)->GetZ();

  }

  fCalcDensity->plasmaE = 28.816 * sqrt(fMaterial->GetDensity()/(g/cm3)

    * sumZ/sumA);

  fCalcDensity->meanexcite = fMeanExcitationEnergy/eV;

  // In case of a fit failure, SetupFermiDeltaCalc prints a warning.

  // We then delete fCalcDensity so it won't be used.

  if(!SetupFermiDeltaCalc(fCalcDensity)){

    free(fCalcDensity->sternf);

    free(fCalcDensity->levE);

    delete fCalcDensity;

    fCalcDensity = nullptr;

    fCalcDensityFailure = true;

  }

}

G4double G4IonisParamMat::DensityCorrection(G4double x)

{

  // If we are set up to calculate the density effect precisely,

  // do that.  Otherwise, fall back to Sternheimer 3-part approximations.

  if(fCalcDensity == nullptr && !fCalcDensityFailure)

    SetSternheimerExactDensityEffect();

  // x = log10(beta*gamma)

  G4double approx = 0.0;

  if(x < fX0density) {

    if(fD0density > 0.0)

      approx = fD0density*G4Exp(twoln10*(x - fX0density));

  }

  else if(x >= fX1density) {

    approx = twoln10*x - fCdensity;

  }

  else{

    approx = twoln10*x - fCdensity

        + fAdensity*G4Exp(G4Log(fX1density - x)*fMdensity);

  }

  if(fCalcDensity != nullptr){

    const double exact = DoFermiDeltaCalc(fCalcDensity, x);

    if(approx > 0 && exact < 0){

      G4cerr << "Error: Sternheimer fit failed for "

        << fMaterial->GetName() << ". x = " << x << ": Delta = "

        << exact << " exact " << approx  << " approx" << G4endl;

      return approx;

    }

    // Fall back to approx if exact and approx are very different, under the

    // assumption that this means the exact calculation has gone haywire

    // somehow, with the exception of the case where approx is negative.  I

    // have seen this clearly-wrong result occur for substances with extremely

    // low density (1e-25 g/cc).

    if(approx >= 0 && fabs(exact - approx) > 1){

      G4cerr << "Error: Sternheimer exact, " << exact << ", and approx, "

        << approx << " are too different for "

        << fMaterial->GetName() << ", x = " << x << G4endl;

      return approx;

    }

    return exact;

  }

  return approx;

}