G4PenelopeComptonModel.cc

Go to the documentation of this file.

//
// ********************************************************************
// * 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.          *
// ********************************************************************
//
//
// Author: Luciano Pandola
//
// History:
// --------
// 15 Feb 2010   L Pandola  Implementation
// 18 Mar 2010   L Pandola  Removed GetAtomsPerMolecule(), now demanded
//                            to G4PenelopeOscillatorManager
// 01 Feb 2011   L Pandola  Suppress fake energy-violation warning when Auger is
//                            active.
//                          Make sure that fluorescence/Auger is generated only
//                            if above threshold
// 24 May 2011   L Pandola  Renamed (make v2008 as default Penelope)
// 10 Jun 2011   L Pandola  Migrate atomic deexcitation interface
// 09 Oct 2013   L Pandola  Migration to MT
// 25 Jul 2023   D Iuso     Fix for possible infinite loops due to FP
//
#include "G4PenelopeComptonModel.hh"
#include "G4PhysicalConstants.hh"
#include "G4SystemOfUnits.hh"
#include "G4ParticleDefinition.hh"
#include "G4MaterialCutsCouple.hh"
#include "G4DynamicParticle.hh"
#include "G4VEMDataSet.hh"
#include "G4PhysicsTable.hh"
#include "G4PhysicsLogVector.hh"
#include "G4AtomicTransitionManager.hh"
#include "G4AtomicShell.hh"
#include "G4Gamma.hh"
#include "G4Electron.hh"
#include "G4PenelopeOscillatorManager.hh"
#include "G4PenelopeOscillator.hh"
#include "G4LossTableManager.hh"
#include "G4Exp.hh"
 
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 
G4PenelopeComptonModel::G4PenelopeComptonModel(const G4ParticleDefinition* part,
                           const G4String& nam)
  :G4VEmModel(nam),fParticleChange(nullptr),fParticle(nullptr),
   fAtomDeexcitation(nullptr),
   fOscManager(nullptr),fIsInitialised(false)
{
  fIntrinsicLowEnergyLimit = 100.0*eV;
  fIntrinsicHighEnergyLimit = 100.0*GeV;
  SetHighEnergyLimit(fIntrinsicHighEnergyLimit);
  //
  fOscManager = G4PenelopeOscillatorManager::GetOscillatorManager();
 
  if (part)
    SetParticle(part);
 
  fVerboseLevel= 0;
  // Verbosity scale:
  // 0 = nothing
  // 1 = warning for energy non-conservation
  // 2 = details of energy budget
  // 3 = calculation of cross sections, file openings, sampling of atoms
  // 4 = entering in methods
 
  //Mark this model as "applicable" for atomic deexcitation
  SetDeexcitationFlag(true);
 
  fTransitionManager = G4AtomicTransitionManager::Instance();
}
 
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 
G4PenelopeComptonModel::~G4PenelopeComptonModel()
{;}
 
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 
void G4PenelopeComptonModel::Initialise(const G4ParticleDefinition* part,
                      const G4DataVector&)
{
  if (fVerboseLevel > 3)
    G4cout << "Calling G4PenelopeComptonModel::Initialise()" << G4endl;
 
  fAtomDeexcitation = G4LossTableManager::Instance()->AtomDeexcitation();
  //Issue warning if the AtomicDeexcitation has not been declared
  if (!fAtomDeexcitation)
    {
      G4cout << G4endl;
      G4cout << "WARNING from G4PenelopeComptonModel " << G4endl;
      G4cout << "Atomic de-excitation module is not instantiated, so there will not be ";
      G4cout << "any fluorescence/Auger emission." << G4endl;
      G4cout << "Please make sure this is intended" << G4endl;
    }
 
  SetParticle(part);
 
  if (IsMaster() && part == fParticle)
    {
 
      if (fVerboseLevel > 0)
    {
      G4cout << "Penelope Compton model v2008 is initialized " << G4endl
         << "Energy range: "
         << LowEnergyLimit() / keV << " keV - "
         << HighEnergyLimit() / GeV << " GeV";
    }
      //Issue a warning, if the model is going to be used down to a
      //energy which is outside the validity of the model itself
      if (LowEnergyLimit() < fIntrinsicLowEnergyLimit)
    {
      G4ExceptionDescription ed;
      ed << "Using the Penelope Compton model outside its intrinsic validity range. "
         << G4endl;
      ed << "-> LowEnergyLimit() in process = " << LowEnergyLimit()/keV << "keV " << G4endl;
      ed << "-> Instrinsic low-energy limit = " << fIntrinsicLowEnergyLimit/keV << "keV "
         << G4endl;
      ed << "Result of the simulation have to be taken with care" << G4endl;
      G4Exception("G4PenelopeComptonModel::Initialise()",
              "em2100",JustWarning,ed);
    }
    }
 
  if(fIsInitialised) return;
  fParticleChange = GetParticleChangeForGamma();
  fIsInitialised = true;
 
}
 
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 
void G4PenelopeComptonModel::InitialiseLocal(const G4ParticleDefinition* part,
                                                     G4VEmModel *masterModel)
{
  if (fVerboseLevel > 3)
    G4cout << "Calling  G4PenelopeComptonModel::InitialiseLocal()" << G4endl;
  //
  //Check that particle matches: one might have multiple master models (e.g.
  //for e+ and e-).
  //
  if (part == fParticle)
    {
      //Get the const table pointers from the master to the workers
      const G4PenelopeComptonModel* theModel =
        static_cast<G4PenelopeComptonModel*> (masterModel);
 
      //Same verbosity for all workers, as the master
      fVerboseLevel = theModel->fVerboseLevel;
    }
  return;
}
 
 
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 
G4double G4PenelopeComptonModel::CrossSectionPerVolume(const G4Material* material,
                                               const G4ParticleDefinition* p,
                                               G4double energy,
                                               G4double,
                                               G4double)
{
  // Penelope model v2008 to calculate the Compton scattering cross section:
  // D. Brusa et al., Nucl. Instrum. Meth. A 379 (1996) 167
  //
  // The cross section for Compton scattering is calculated according to the Klein-Nishina
  // formula for energy > 5 MeV.
  // For E < 5 MeV it is used a parametrization for the differential cross-section dSigma/dOmega,
  // which is integrated numerically in cos(theta), G4PenelopeComptonModel::DifferentialCrossSection().
  // The parametrization includes the J(p)
  // distribution profiles for the atomic shells, that are tabulated from Hartree-Fock calculations
  // from F. Biggs et al., At. Data Nucl. Data Tables 16 (1975) 201
  //
  if (fVerboseLevel > 3)
    G4cout << "Calling CrossSectionPerVolume() of G4PenelopeComptonModel" << G4endl;
  SetupForMaterial(p, material, energy);
 
  G4double cs = 0;
  //Force null cross-section if below the low-energy edge of the table
  if (energy < LowEnergyLimit())
    return cs;
 
  //Retrieve the oscillator table for this material
  G4PenelopeOscillatorTable* theTable = fOscManager->GetOscillatorTableCompton(material);
 
  if (energy < 5*MeV) //explicit calculation for E < 5 MeV
    {
      size_t numberOfOscillators = theTable->size();
      for (size_t i=0;i<numberOfOscillators;i++)
    {
      G4PenelopeOscillator* theOsc = (*theTable)[i];
      //sum contributions coming from each oscillator
      cs += OscillatorTotalCrossSection(energy,theOsc);
    }
    }
  else //use Klein-Nishina for E>5 MeV
    cs = KleinNishinaCrossSection(energy,material);
 
  //cross sections are in units of pi*classic_electr_radius^2
  cs *= pi*classic_electr_radius*classic_electr_radius;
 
  //Now, cs is the cross section *per molecule*, let's calculate the
  //cross section per volume
  G4double atomDensity = material->GetTotNbOfAtomsPerVolume();
  G4double atPerMol =  fOscManager->GetAtomsPerMolecule(material);
 
  if (fVerboseLevel > 3)
    G4cout << "Material " << material->GetName() << " has " << atPerMol <<
      "atoms per molecule" << G4endl;
 
  G4double moleculeDensity = 0.;
 
  if (atPerMol)
    moleculeDensity = atomDensity/atPerMol;
 
  G4double csvolume = cs*moleculeDensity;
 
  if (fVerboseLevel > 2)
    G4cout << "Compton mean free path at " << energy/keV << " keV for material " <<
            material->GetName() << " = " << (1./csvolume)/mm << " mm" << G4endl;
  return csvolume;
}
 
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 
//This is a dummy method. Never inkoved by the tracking, it just issues
//a warning if one tries to get Cross Sections per Atom via the
//G4EmCalculator.
G4double G4PenelopeComptonModel::ComputeCrossSectionPerAtom(const G4ParticleDefinition*,
                                                             G4double,
                                                             G4double,
                                                             G4double,
                                                             G4double,
                                                             G4double)
{
  G4cout << "*** G4PenelopeComptonModel -- WARNING ***" << G4endl;
  G4cout << "Penelope Compton model v2008 does not calculate cross section _per atom_ " << G4endl;
  G4cout << "so the result is always zero. For physics values, please invoke " << G4endl;
  G4cout << "GetCrossSectionPerVolume() or GetMeanFreePath() via the G4EmCalculator" << G4endl;
  return 0;
}
 
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 
void G4PenelopeComptonModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
                           const G4MaterialCutsCouple* couple,
                          const G4DynamicParticle* aDynamicGamma,
                          G4double,
                          G4double)
{
  // Penelope model v2008 to sample the Compton scattering final state.
  // D. Brusa et al., Nucl. Instrum. Meth. A 379 (1996) 167
  // The model determines also the original shell from which the electron is expelled,
  // in order to produce fluorescence de-excitation (from G4DeexcitationManager)
  //
  // The final state for Compton scattering is calculated according to the Klein-Nishina
  // formula for energy > 5 MeV. In this case, the Doppler broadening is negligible and
  // one can assume that the target electron is at rest.
  // For E < 5 MeV it is used the parametrization for the differential cross-section dSigma/dOmega,
  // to sample the scattering angle and the energy of the emerging electron, which is
  // G4PenelopeComptonModel::DifferentialCrossSection(). The rejection method is
  // used to sample cos(theta). The efficiency increases monotonically with photon energy and is
  // nearly independent on the Z; typical values are 35%, 80% and 95% for 1 keV, 1 MeV and 10 MeV,
  // respectively.
  // The parametrization includes the J(p) distribution profiles for the atomic shells, that are
  // tabulated
  // from Hartree-Fock calculations from F. Biggs et al., At. Data Nucl. Data Tables 16 (1975) 201.
  // Doppler broadening is included.
  //
 
  if (fVerboseLevel > 3)
    G4cout << "Calling SampleSecondaries() of G4PenelopeComptonModel" << G4endl;
 
  G4double photonEnergy0 = aDynamicGamma->GetKineticEnergy();
 
  // do nothing below the threshold
  // should never get here because the XS is zero below the limit
  if(photonEnergy0 < LowEnergyLimit())
    return;
 
  G4ParticleMomentum photonDirection0 = aDynamicGamma->GetMomentumDirection();
  const G4Material* material = couple->GetMaterial();
 
  G4PenelopeOscillatorTable* theTable = fOscManager->GetOscillatorTableCompton(material);
 
  const G4int nmax = 64;
  G4double rn[nmax]={0.0};
  G4double pac[nmax]={0.0};
 
  G4double S=0.0;
  G4double epsilon = 0.0;
  G4double cosTheta = 1.0;
  G4double hartreeFunc = 0.0;
  G4double oscStren = 0.0;
  size_t numberOfOscillators = theTable->size();
  size_t targetOscillator = 0;
  G4double ionEnergy = 0.0*eV;
 
  G4double ek = photonEnergy0/electron_mass_c2;
  G4double ek2 = 2.*ek+1.0;
  G4double eks = ek*ek;
  G4double ek1 = eks-ek2-1.0;
 
  G4double taumin = 1.0/ek2;
  //This is meant to fix a possible (rare) floating point exception in the sampling of tau below,
  //causing an infinite loop. The maximum of tau is not 1., but the closest double which can 
  //be represented (i.e. ~ 1. - 1e-16). Fix by Domenico Iuso 
  static G4double taumax = std::nexttoward(1.0,0.0);
  if (fVerboseLevel > 3)
    G4cout << "G4PenelopeComptonModel: maximum value of tau: 1 - " << 1.-taumax << G4endl;
  //To here. 
  G4double a1 = G4Log(ek2);
  G4double a2 = a1+2.0*ek*(1.0+ek)/(ek2*ek2);
 
  G4double TST = 0;
  G4double tau = 0.;
 
  //If the incoming photon is above 5 MeV, the quicker approach based on the
  //pure Klein-Nishina formula is used
  if (photonEnergy0 > 5*MeV)
    {
      do{
    do{
      if ((a2*G4UniformRand()) < a1)
        tau = std::pow(taumin,G4UniformRand());
      else
        tau = std::sqrt(1.0+G4UniformRand()*(taumin*taumin-1.0));
      //rejection function
      TST = (1.0+tau*(ek1+tau*(ek2+tau*eks)))/(eks*tau*(1.0+tau*tau));
    }while (G4UniformRand()> TST);
    if (tau > taumax) tau = taumax; //prevent FP exception causing infinite loop
    epsilon=tau;
    cosTheta = 1.0 - (1.0-tau)/(ek*tau);
 
    //Target shell electrons
    TST = fOscManager->GetTotalZ(material)*G4UniformRand();
    targetOscillator = numberOfOscillators-1; //last level
    S=0.0;
    G4bool levelFound = false;
    for (size_t j=0;j<numberOfOscillators && !levelFound; j++)
      {
        S += (*theTable)[j]->GetOscillatorStrength();
        if (S > TST)
          {
        targetOscillator = j;
        levelFound = true;
          }
      }
    //check whether the level is valid
    ionEnergy = (*theTable)[targetOscillator]->GetIonisationEnergy();
      }while((epsilon*photonEnergy0-photonEnergy0+ionEnergy) >0);
    }
  else //photonEnergy0 < 5 MeV
    {
      //Incoherent scattering function for theta=PI
      G4double s0=0.0;
      G4double pzomc=0.0;
      G4double rni=0.0;
      G4double aux=0.0;
      for (size_t i=0;i<numberOfOscillators;i++)
    {
      ionEnergy = (*theTable)[i]->GetIonisationEnergy();
      if (photonEnergy0 > ionEnergy)
        {
          G4double aux2 = photonEnergy0*(photonEnergy0-ionEnergy)*2.0;
          hartreeFunc = (*theTable)[i]->GetHartreeFactor();
          oscStren = (*theTable)[i]->GetOscillatorStrength();
          pzomc = hartreeFunc*(aux2-electron_mass_c2*ionEnergy)/
        (electron_mass_c2*std::sqrt(2.0*aux2+ionEnergy*ionEnergy));
          if (pzomc > 0)
        rni = 1.0-0.5*G4Exp(0.5-(std::sqrt(0.5)+std::sqrt(2.0)*pzomc)*
                       (std::sqrt(0.5)+std::sqrt(2.0)*pzomc));
          else
        rni = 0.5*G4Exp(0.5-(std::sqrt(0.5)-std::sqrt(2.0)*pzomc)*
                   (std::sqrt(0.5)-std::sqrt(2.0)*pzomc));
          s0 += oscStren*rni;
        }
    }
      //Sampling tau
      G4double cdt1 = 0.;
      do
    {
      if ((G4UniformRand()*a2) < a1)
        tau = std::pow(taumin,G4UniformRand());
      else
        tau = std::sqrt(1.0+G4UniformRand()*(taumin*taumin-1.0));
      if (tau > taumax) tau = taumax; //prevent FP exception causing infinite loop
      cdt1 = (1.0-tau)/(ek*tau);
      //Incoherent scattering function
      S = 0.;
      for (size_t i=0;i<numberOfOscillators;i++)
        {
          ionEnergy = (*theTable)[i]->GetIonisationEnergy();
          if (photonEnergy0 > ionEnergy) //sum only on excitable levels
        {
          aux = photonEnergy0*(photonEnergy0-ionEnergy)*cdt1;
          hartreeFunc = (*theTable)[i]->GetHartreeFactor();
          oscStren = (*theTable)[i]->GetOscillatorStrength();
          pzomc = hartreeFunc*(aux-electron_mass_c2*ionEnergy)/
            (electron_mass_c2*std::sqrt(2.0*aux+ionEnergy*ionEnergy));
          if (pzomc > 0)
            rn[i] = 1.0-0.5*G4Exp(0.5-(std::sqrt(0.5)+std::sqrt(2.0)*pzomc)*
                         (std::sqrt(0.5)+std::sqrt(2.0)*pzomc));
          else
            rn[i] = 0.5*G4Exp(0.5-(std::sqrt(0.5)-std::sqrt(2.0)*pzomc)*
                     (std::sqrt(0.5)-std::sqrt(2.0)*pzomc));
          S += oscStren*rn[i];
          pac[i] = S;
        }
          else
        pac[i] = S-1e-6;
        }
      //Rejection function
      TST = S*(1.0+tau*(ek1+tau*(ek2+tau*eks)))/(eks*tau*(1.0+tau*tau));
    }while ((G4UniformRand()*s0) > TST);
 
      cosTheta = 1.0 - cdt1;
      G4double fpzmax=0.0,fpz=0.0;
      G4double A=0.0;
      //Target electron shell
      do
    {
      do
        {
          TST = S*G4UniformRand();
          targetOscillator = numberOfOscillators-1; //last level
          G4bool levelFound = false;
          for (size_t i=0;i<numberOfOscillators && !levelFound;i++)
        {
          if (pac[i]>TST)
            {
              targetOscillator = i;
              levelFound = true;
            }
        }
          A = G4UniformRand()*rn[targetOscillator];
          hartreeFunc = (*theTable)[targetOscillator]->GetHartreeFactor();
          oscStren = (*theTable)[targetOscillator]->GetOscillatorStrength();
          if (A < 0.5)
        pzomc = (std::sqrt(0.5)-std::sqrt(0.5-G4Log(2.0*A)))/
          (std::sqrt(2.0)*hartreeFunc);
          else
        pzomc = (std::sqrt(0.5-G4Log(2.0-2.0*A))-std::sqrt(0.5))/
          (std::sqrt(2.0)*hartreeFunc);
        } while (pzomc < -1);
 
      // F(EP) rejection
      G4double XQC = 1.0+tau*(tau-2.0*cosTheta);
      G4double AF = std::sqrt(XQC)*(1.0+tau*(tau-cosTheta)/XQC);
      if (AF > 0)
        fpzmax = 1.0+AF*0.2;
      else
        fpzmax = 1.0-AF*0.2;
      fpz = 1.0+AF*std::max(std::min(pzomc,0.2),-0.2);
    }while ((fpzmax*G4UniformRand())>fpz);
 
      //Energy of the scattered photon
      G4double T = pzomc*pzomc;
      G4double b1 = 1.0-T*tau*tau;
      G4double b2 = 1.0-T*tau*cosTheta;
      if (pzomc > 0.0)
    epsilon = (tau/b1)*(b2+std::sqrt(std::abs(b2*b2-b1*(1.0-T))));
      else
    epsilon = (tau/b1)*(b2-std::sqrt(std::abs(b2*b2-b1*(1.0-T))));
    } //energy < 5 MeV
 
  //Ok, the kinematics has been calculated.
  G4double sinTheta = std::sqrt(1-cosTheta*cosTheta);
  G4double phi = twopi * G4UniformRand() ;
  G4double dirx = sinTheta * std::cos(phi);
  G4double diry = sinTheta * std::sin(phi);
  G4double dirz = cosTheta ;
 
  // Update G4VParticleChange for the scattered photon
  G4ThreeVector photonDirection1(dirx,diry,dirz);
  photonDirection1.rotateUz(photonDirection0);
  fParticleChange->ProposeMomentumDirection(photonDirection1) ;
 
  G4double photonEnergy1 = epsilon * photonEnergy0;
 
  if (photonEnergy1 > 0.)
    fParticleChange->SetProposedKineticEnergy(photonEnergy1) ;
  else
  {
    fParticleChange->SetProposedKineticEnergy(0.) ;
    fParticleChange->ProposeTrackStatus(fStopAndKill);
  }
 
  //Create scattered electron
  G4double diffEnergy = photonEnergy0*(1-epsilon);
  ionEnergy = (*theTable)[targetOscillator]->GetIonisationEnergy();
 
  G4double Q2 =
    photonEnergy0*photonEnergy0+photonEnergy1*(photonEnergy1-2.0*photonEnergy0*cosTheta);
  G4double cosThetaE = 0.; //scattering angle for the electron
 
  if (Q2 > 1.0e-12)
    cosThetaE = (photonEnergy0-photonEnergy1*cosTheta)/std::sqrt(Q2);
  else
    cosThetaE = 1.0;
  G4double sinThetaE = std::sqrt(1-cosThetaE*cosThetaE);
 
  //Now, try to handle fluorescence
  //Notice: merged levels are indicated with Z=0 and flag=30
  G4int shFlag = (*theTable)[targetOscillator]->GetShellFlag();
  G4int Z = (G4int) (*theTable)[targetOscillator]->GetParentZ();
 
  //initialize here, then check photons created by Atomic-Deexcitation, and the final state e-
  G4double bindingEnergy = 0.*eV;
  const G4AtomicShell* shell = 0;
 
  //Real level
  if (Z > 0 && shFlag<30)
    {
      shell = fTransitionManager->Shell(Z,shFlag-1);
      bindingEnergy = shell->BindingEnergy();
    }
 
  G4double ionEnergyInPenelopeDatabase = ionEnergy;
  //protection against energy non-conservation
  ionEnergy = std::max(bindingEnergy,ionEnergyInPenelopeDatabase);
 
  //subtract the excitation energy. If not emitted by fluorescence
  //the ionization energy is deposited as local energy deposition
  G4double eKineticEnergy = diffEnergy - ionEnergy;
  G4double localEnergyDeposit = ionEnergy;
  G4double energyInFluorescence = 0.; //testing purposes only
  G4double energyInAuger = 0; //testing purposes
 
  if (eKineticEnergy < 0)
    {
      //It means that there was some problem/mismatch between the two databases.
      //Try to make it work
      //In this case available Energy (diffEnergy) < ionEnergy
      //Full residual energy is deposited locally
      localEnergyDeposit = diffEnergy;
      eKineticEnergy = 0.0;
    }
 
  //the local energy deposit is what remains: part of this may be spent for fluorescence.
  //Notice: shell might be nullptr (invalid!) if shFlag=30. Must be protected
  //Now, take care of fluorescence, if required
  if (fAtomDeexcitation && shell)
    {
      G4int index = couple->GetIndex();
      if (fAtomDeexcitation->CheckDeexcitationActiveRegion(index))
    {
      size_t nBefore = fvect->size();
      fAtomDeexcitation->GenerateParticles(fvect,shell,Z,index);
      size_t nAfter = fvect->size();
 
      if (nAfter > nBefore) //actual production of fluorescence
        {
          for (size_t j=nBefore;j<nAfter;j++) //loop on products
        {
          G4double itsEnergy = ((*fvect)[j])->GetKineticEnergy();
          if (itsEnergy < localEnergyDeposit) // valid secondary, generate it
            {
              localEnergyDeposit -= itsEnergy;
              if (((*fvect)[j])->GetParticleDefinition() == G4Gamma::Definition())
            energyInFluorescence += itsEnergy;
              else if (((*fvect)[j])->GetParticleDefinition() == 
                   G4Electron::Definition())
            energyInAuger += itsEnergy;
            }
          else //invalid secondary: takes more than the available energy: delete it
            {
              delete (*fvect)[j];
              (*fvect)[j] = nullptr;
            }
        }
        }
      
    }
    }
 
  //Always produce explicitly the electron
  G4DynamicParticle* electron = 0;
 
  G4double xEl = sinThetaE * std::cos(phi+pi);
  G4double yEl = sinThetaE * std::sin(phi+pi);
  G4double zEl = cosThetaE;
  G4ThreeVector eDirection(xEl,yEl,zEl); //electron direction
  eDirection.rotateUz(photonDirection0);
  electron = new G4DynamicParticle (G4Electron::Electron(),
                    eDirection,eKineticEnergy) ;
  fvect->push_back(electron);
 
  if (localEnergyDeposit < 0) //Should not be: issue a G4Exception (warning)
    {
      G4Exception("G4PenelopeComptonModel::SampleSecondaries()",
           "em2099",JustWarning,"WARNING: Negative local energy deposit");
      localEnergyDeposit=0.;
    }
  fParticleChange->ProposeLocalEnergyDeposit(localEnergyDeposit);
 
  G4double electronEnergy = 0.;
  if (electron)
    electronEnergy = eKineticEnergy;
  if (fVerboseLevel > 1)
    {
      G4cout << "-----------------------------------------------------------" << G4endl;
      G4cout << "Energy balance from G4PenelopeCompton" << G4endl;
      G4cout << "Incoming photon energy: " << photonEnergy0/keV << " keV" << G4endl;
      G4cout << "-----------------------------------------------------------" << G4endl;
      G4cout << "Scattered photon: " << photonEnergy1/keV << " keV" << G4endl;
      G4cout << "Scattered electron " << electronEnergy/keV << " keV" << G4endl;
      if (energyInFluorescence)
    G4cout << "Fluorescence x-rays: " << energyInFluorescence/keV << " keV" << G4endl;
      if (energyInAuger)
    G4cout << "Auger electrons: " << energyInAuger/keV << " keV" << G4endl;
      G4cout << "Local energy deposit " << localEnergyDeposit/keV << " keV" << G4endl;
      G4cout << "Total final state: " << (photonEnergy1+electronEnergy+energyInFluorescence+
                      localEnergyDeposit+energyInAuger)/keV <<
    " keV" << G4endl;
      G4cout << "-----------------------------------------------------------" << G4endl;
    }
  if (fVerboseLevel > 0)
    {
      G4double energyDiff = std::fabs(photonEnergy1+
                      electronEnergy+energyInFluorescence+
                      localEnergyDeposit+energyInAuger-photonEnergy0);
      if (energyDiff > 0.05*keV)
    G4cout << "Warning from G4PenelopeCompton: problem with energy conservation: " <<
      (photonEnergy1+electronEnergy+energyInFluorescence+energyInAuger+localEnergyDeposit)/keV <<
      " keV (final) vs. " <<
      photonEnergy0/keV << " keV (initial)" << G4endl;
    }
}
 
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 
G4double G4PenelopeComptonModel::DifferentialCrossSection(G4double cosTheta,G4double energy,
                                G4PenelopeOscillator* osc)
{
  //
  // Penelope model v2008. Single differential cross section *per electron*
  // for photon Compton scattering by
  // electrons in the given atomic oscillator, differential in the direction of the
  // scattering photon. This is in units of pi*classic_electr_radius**2
  //
  // D. Brusa et al., Nucl. Instrum. Meth. A 379 (1996) 167
  // The parametrization includes the J(p) distribution profiles for the atomic shells,
  // that are tabulated from Hartree-Fock calculations
  // from F. Biggs et al., At. Data Nucl. Data Tables 16 (1975) 201
  //
  G4double ionEnergy = osc->GetIonisationEnergy();
  G4double harFunc = osc->GetHartreeFactor();
 
  static const G4double k2 = std::sqrt(2.);
  static const G4double k1 = 1./k2;
 
  if (energy < ionEnergy)
    return 0;
 
  //energy of the Compton line
  G4double cdt1 = 1.0-cosTheta;
  G4double EOEC = 1.0+(energy/electron_mass_c2)*cdt1;
  G4double ECOE = 1.0/EOEC;
 
  //Incoherent scattering function (analytical profile)
  G4double aux = energy*(energy-ionEnergy)*cdt1;
  G4double Pzimax =
    (aux - electron_mass_c2*ionEnergy)/(electron_mass_c2*std::sqrt(2*aux+ionEnergy*ionEnergy));
  G4double sia = 0.0;
  G4double x = harFunc*Pzimax;
  if (x > 0)
    sia = 1.0-0.5*G4Exp(0.5-(k1+k2*x)*(k1+k2*x));
  else
    sia = 0.5*G4Exp(0.5-(k1-k2*x)*(k1-k2*x));
 
  //1st order correction, integral of Pz times the Compton profile.
  //Calculated approximately using a free-electron gas profile
  G4double pf = 3.0/(4.0*harFunc);
  if (std::fabs(Pzimax) < pf)
    {
      G4double QCOE2 = 1.0+ECOE*ECOE-2.0*ECOE*cosTheta;
      G4double p2 = Pzimax*Pzimax;
      G4double dspz = std::sqrt(QCOE2)*
    (1.0+ECOE*(ECOE-cosTheta)/QCOE2)*harFunc
    *0.25*(2*p2-(p2*p2)/(pf*pf)-(pf*pf));
      sia += std::max(dspz,-1.0*sia);
    }
 
  G4double XKN = EOEC+ECOE-1.0+cosTheta*cosTheta;
 
  //Differential cross section (per electron, in units of pi*classic_electr_radius**2)
  G4double diffCS = ECOE*ECOE*XKN*sia;
 
  return diffCS;
}
 
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 
G4double G4PenelopeComptonModel::OscillatorTotalCrossSection(G4double energy,G4PenelopeOscillator* osc)
{
  //Total cross section (integrated) for the given oscillator in units of
  //pi*classic_electr_radius^2
 
  //Integrate differential cross section for each oscillator
  G4double stre = osc->GetOscillatorStrength();
 
  // here one uses the  using the 20-point
  // Gauss quadrature method with an adaptive bipartition scheme
  const G4int npoints=10;
  const G4int ncallsmax=20000;
  const G4int nst=256;
  static const G4double Abscissas[10] = {7.652651133497334e-02,2.2778585114164508e-01,3.7370608871541956e-01,
                  5.1086700195082710e-01,6.3605368072651503e-01,7.4633190646015079e-01,
                  8.3911697182221882e-01,9.1223442825132591e-01,9.6397192727791379e-01,
                  9.9312859918509492e-01};
  static const G4double Weights[10] = {1.5275338713072585e-01,1.4917298647260375e-01,1.4209610931838205e-01,
                1.3168863844917663e-01,1.1819453196151842e-01,1.0193011981724044e-01,
                8.3276741576704749e-02,6.2672048334109064e-02,4.0601429800386941e-02,
                1.7614007139152118e-02};
 
  G4double MaxError = 1e-5;
  //Error control
  G4double Ctol = std::min(std::max(MaxError,1e-13),1e-02);
  G4double Ptol = 0.01*Ctol;
  G4double Err=1e35;
 
  //Gauss integration from -1 to 1
  G4double LowPoint = -1.0;
  G4double HighPoint = 1.0;
 
  G4double h=HighPoint-LowPoint;
  G4double sumga=0.0;
  G4double a=0.5*(HighPoint-LowPoint);
  G4double b=0.5*(HighPoint+LowPoint);
  G4double c=a*Abscissas[0];
  G4double d= Weights[0]*
    (DifferentialCrossSection(b+c,energy,osc)+DifferentialCrossSection(b-c,energy,osc));
  for (G4int i=2;i<=npoints;i++)
    {
      c=a*Abscissas[i-1];
      d += Weights[i-1]*
    (DifferentialCrossSection(b+c,energy,osc)+DifferentialCrossSection(b-c,energy,osc));
    }
  G4int icall = 2*npoints;
  G4int LH=1;
  G4double S[nst],x[nst],sn[nst],xrn[nst];
  S[0]=d*a;
  x[0]=LowPoint;
 
  G4bool loopAgain = true;
 
  //Adaptive bipartition scheme
  do{
    G4double h0=h;
    h=0.5*h; //bipartition
    G4double sumr=0;
    G4int LHN=0;
    G4double si,xa,xb,xc;
    for (G4int i=1;i<=LH;i++){
      si=S[i-1];
      xa=x[i-1];
      xb=xa+h;
      xc=xa+h0;
      a=0.5*(xb-xa);
      b=0.5*(xb+xa);
      c=a*Abscissas[0];
      G4double dLocal = Weights[0]*
    (DifferentialCrossSection(b+c,energy,osc)+DifferentialCrossSection(b-c,energy,osc));
 
      for (G4int j=1;j<npoints;j++)
    {
      c=a*Abscissas[j];
      dLocal += Weights[j]*
        (DifferentialCrossSection(b+c,energy,osc)+DifferentialCrossSection(b-c,energy,osc));
    }
      G4double s1=dLocal*a;
      a=0.5*(xc-xb);
      b=0.5*(xc+xb);
      c=a*Abscissas[0];
      dLocal=Weights[0]*
    (DifferentialCrossSection(b+c,energy,osc)+DifferentialCrossSection(b-c,energy,osc));
 
      for (G4int j=1;j<npoints;j++)
    {
      c=a*Abscissas[j];
      dLocal += Weights[j]*
        (DifferentialCrossSection(b+c,energy,osc)+DifferentialCrossSection(b-c,energy,osc));
    }
      G4double s2=dLocal*a;
      icall=icall+4*npoints;
      G4double s12=s1+s2;
      if (std::abs(s12-si)<std::max(Ptol*std::abs(s12),1e-35))
    sumga += s12;
      else
    {
      sumr += s12;
      LHN += 2;
      sn[LHN-1]=s2;
      xrn[LHN-1]=xb;
      sn[LHN-2]=s1;
      xrn[LHN-2]=xa;
    }
 
      if (icall>ncallsmax || LHN>nst)
    {
      G4cout << "G4PenelopeComptonModel: " << G4endl;
      G4cout << "LowPoint: " << LowPoint << ", High Point: " << HighPoint << G4endl;
      G4cout << "Tolerance: " << MaxError << G4endl;
      G4cout << "Calls: " << icall << ", Integral: " << sumga << ", Error: " << Err << G4endl;
      G4cout << "Number of open subintervals: " << LHN << G4endl;
      G4cout << "WARNING: the required accuracy has not been attained" << G4endl;
      loopAgain = false;
    }
    }
    Err=std::abs(sumr)/std::max(std::abs(sumr+sumga),1e-35);
    if (Err < Ctol || LHN == 0)
      loopAgain = false; //end of cycle
    LH=LHN;
    for (G4int i=0;i<LH;i++)
      {
    S[i]=sn[i];
    x[i]=xrn[i];
      }
  }while(Ctol < 1.0 && loopAgain);
 
  G4double xs = stre*sumga;
 
  return xs;
}
 
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 
G4double G4PenelopeComptonModel::KleinNishinaCrossSection(G4double energy,
                           const G4Material* material)
{
  // use Klein-Nishina formula
  // total cross section in units of pi*classic_electr_radius^2
  G4double cs = 0;
 
  G4double ek =energy/electron_mass_c2;
  G4double eks = ek*ek;
  G4double ek2 = 1.0+ek+ek;
  G4double ek1 = eks-ek2-1.0;
 
  G4double t0 = 1.0/ek2;
  G4double csl = 0.5*eks*t0*t0+ek2*t0+ek1*G4Log(t0)-(1.0/t0);
 
  G4PenelopeOscillatorTable* theTable = fOscManager->GetOscillatorTableCompton(material);
 
  for (size_t i=0;i<theTable->size();i++)
    {
      G4PenelopeOscillator* theOsc = (*theTable)[i];
      G4double ionEnergy = theOsc->GetIonisationEnergy();
      G4double tau=(energy-ionEnergy)/energy;
      if (tau > t0)
    {
      G4double csu = 0.5*eks*tau*tau+ek2*tau+ek1*G4Log(tau)-(1.0/tau);
      G4double stre = theOsc->GetOscillatorStrength();
 
      cs += stre*(csu-csl);
    }
    }
  cs /= (ek*eks);
 
  return cs;
 
}
 
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...
 
void G4PenelopeComptonModel::SetParticle(const G4ParticleDefinition* p)
{
  if(!fParticle) {
    fParticle = p;
  }
}