//
// ********************************************************************
// * 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.          *
// ********************************************************************
//
// $Id: G4Penelope08ComptonModel.cc,v 1.6 2010/04/12 13:53:29 pandola Exp $
// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
//
// Author: Luciano Pandola
//
// History:
// --------
// 15 Feb 2010   L Pandola    Implementation
// 18 Mar 2010   L. Pandola   Removed GetAtomsPerMolecule(), now demanded 
//                            to G4PenelopeOscillatorManager
//
#include "G4Penelope08ComptonModel.hh"
#include "G4ParticleDefinition.hh"
#include "G4MaterialCutsCouple.hh"
#include "G4ProductionCutsTable.hh"
#include "G4DynamicParticle.hh"
#include "G4VEMDataSet.hh"
#include "G4PhysicsTable.hh"
#include "G4ElementTable.hh"
#include "G4Element.hh"
#include "G4PhysicsLogVector.hh"
#include "G4AtomicTransitionManager.hh"
#include "G4AtomicShell.hh"
#include "G4Gamma.hh"
#include "G4Electron.hh"
#include "G4PenelopeOscillatorManager.hh"
#include "G4PenelopeOscillator.hh"

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....


G4Penelope08ComptonModel::G4Penelope08ComptonModel(const G4ParticleDefinition*,
                                             const G4String& nam)
  :G4VEmModel(nam),isInitialised(false),oscManager(0)
{
  fIntrinsicLowEnergyLimit = 100.0*eV;
  fIntrinsicHighEnergyLimit = 100.0*GeV;
  //  SetLowEnergyLimit(fIntrinsicLowEnergyLimit);
  SetHighEnergyLimit(fIntrinsicHighEnergyLimit);
  //
  oscManager = G4PenelopeOscillatorManager::GetOscillatorManager();

  verboseLevel= 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

  //by default, the model will use atomic deexcitation
  SetDeexcitationFlag(true);
  ActivateAuger(false);

}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....

G4Penelope08ComptonModel::~G4Penelope08ComptonModel()
{;}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....

void G4Penelope08ComptonModel::Initialise(const G4ParticleDefinition*,
					  const G4DataVector&)
{
  if (verboseLevel > 3)
    G4cout << "Calling G4Penelope08ComptonModel::Initialise()" << G4endl;

  if (verboseLevel > 0) {
    G4cout << "Penelope Compton model is initialized " << G4endl
	   << "Energy range: "
	   << LowEnergyLimit() / keV << " keV - "
	   << HighEnergyLimit() / GeV << " GeV"
	   << G4endl;
  }

  if(isInitialised) return;
  fParticleChange = GetParticleChangeForGamma();
  isInitialised = true; 
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....

G4double G4Penelope08ComptonModel::CrossSectionPerVolume(const G4Material* material,
                                               const G4ParticleDefinition* p,
                                               G4double energy,
                                               G4double,
                                               G4double)
{
  // Penelope model 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), G4Penelope08ComptonModel::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 (verboseLevel > 3)
    G4cout << "Calling CrossSectionPerVolume() of G4Penelope08ComptonModel" << G4endl;
 SetupForMaterial(p, material, energy);

  //Retrieve the oscillator table for this material
  G4PenelopeOscillatorTable* theTable = oscManager->GetOscillatorTableCompton(material);

  G4double cs = 0;

  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 =  oscManager->GetAtomsPerMolecule(material);

  if (verboseLevel > 3)
    G4cout << "Material " << material->GetName() << " has " << atPerMol << 
      "atoms per molecule" << G4endl;

  G4double moleculeDensity = 0.;

  if (atPerMol)
    moleculeDensity = atomDensity/atPerMol;

  G4double csvolume = cs*moleculeDensity;
  
  if (verboseLevel > 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 G4Penelope08ComptonModel::ComputeCrossSectionPerAtom(const G4ParticleDefinition*,
                                                             G4double,
                                                             G4double,
                                                             G4double,
                                                             G4double,
                                                             G4double)
{
  G4cout << "*** G4Penelope08ComptonModel -- WARNING ***" << G4endl;
  G4cout << "Penelope Compton model 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 G4Penelope08ComptonModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
					      const G4MaterialCutsCouple* couple,
					      const G4DynamicParticle* aDynamicGamma,
					      G4double,
					      G4double)
{
  
  // Penelope model 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  
  // G4Penelope08ComptonModel::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 (verboseLevel > 3)
    G4cout << "Calling SampleSecondaries() of G4Penelope08ComptonModel" << G4endl;
  
  G4double photonEnergy0 = aDynamicGamma->GetKineticEnergy();

  if (photonEnergy0 <= fIntrinsicLowEnergyLimit)
    {
      fParticleChange->ProposeTrackStatus(fStopAndKill);
      fParticleChange->SetProposedKineticEnergy(0.);
      fParticleChange->ProposeLocalEnergyDeposit(photonEnergy0);
      return ;
    }

  G4ParticleMomentum photonDirection0 = aDynamicGamma->GetMomentumDirection();
  const G4Material* material = couple->GetMaterial();

  G4PenelopeOscillatorTable* theTable = oscManager->GetOscillatorTableCompton(material); 

  const G4int nmax = 64;
  G4double rn[nmax],pac[nmax];
  
  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;
  G4double a1 = std::log(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);
	epsilon=tau;
	cosTheta = 1.0 - (1.0-tau)/(ek*tau);

	//Target shell electrons	
	TST = oscManager->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 aux = photonEnergy0*(photonEnergy0-ionEnergy)*2.0;
	      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) 	
		rni = 1.0-0.5*std::exp(0.5-(std::sqrt(0.5)+std::sqrt(2.0)*pzomc)*
				       (std::sqrt(0.5)+std::sqrt(2.0)*pzomc));		
	      else		  
		rni = 0.5*std::exp(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));	    
	  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*std::exp(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*std::exp(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-std::log(2.0*A)))/
		  (std::sqrt(2.0)*hartreeFunc);	      
	      else		
		pzomc = (std::sqrt(0.5-std::log(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-
  std::vector<G4DynamicParticle*>* photonVector=0;  
  const G4AtomicTransitionManager* transitionManager = G4AtomicTransitionManager::Instance();
  G4double bindingEnergy = 0.*eV;
  G4int shellId = 0;

  //Real level
  if (Z > 0 && shFlag<30)
    {
      const G4AtomicShell* shell = transitionManager->Shell(Z,shFlag-1);
      bindingEnergy = shell->BindingEnergy();
      shellId = shell->ShellId();
    }

  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

  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.
  if(DeexcitationFlag() && Z > 5) {

    const G4ProductionCutsTable* theCoupleTable=
      G4ProductionCutsTable::GetProductionCutsTable();

    size_t index = couple->GetIndex();
    G4double cutg = (*(theCoupleTable->GetEnergyCutsVector(0)))[index];
    G4double cute = (*(theCoupleTable->GetEnergyCutsVector(1)))[index];

    // Generation of fluorescence
    // Data in EADL are available only for Z > 5
    // Protection to avoid generating photons in the unphysical case of
    // shell binding energy > photon energy
    if (localEnergyDeposit > cutg || localEnergyDeposit > cute)
      { 
	G4DynamicParticle* aPhoton;
	deexcitationManager.SetCutForSecondaryPhotons(cutg);
	deexcitationManager.SetCutForAugerElectrons(cute);

	photonVector = deexcitationManager.GenerateParticles(Z,shellId);
	if(photonVector) 
	  {
	    size_t nPhotons = photonVector->size();
	    for (size_t k=0; k<nPhotons; k++)
	      {
		aPhoton = (*photonVector)[k];
		if (aPhoton)
		  {
		    G4double itsEnergy = aPhoton->GetKineticEnergy();
		    if (itsEnergy <= localEnergyDeposit)
		      {
			localEnergyDeposit -= itsEnergy;
			if (aPhoton->GetDefinition() == G4Gamma::Gamma()) 
			  energyInFluorescence += itsEnergy;;
			fvect->push_back(aPhoton);		    
		      }
		    else
		      {
			delete aPhoton;
			(*photonVector)[k]=0;
		      }
		  }
	      }
	    delete photonVector;
	  }
      }
  }

  //Always produce explicitely 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)
    {
      G4cout << "WARNING-" 
	     << "G4Penelope08ComptonModel::SampleSecondaries - Negative energy deposit"
	     << G4endl;
      localEnergyDeposit=0.;
    }
  fParticleChange->ProposeLocalEnergyDeposit(localEnergyDeposit);
  
  G4double electronEnergy = 0.;
  if (verboseLevel > 1)
    {
      G4cout << "-----------------------------------------------------------" << G4endl;
      G4cout << "Energy balance from G4Penelope08Compton" << G4endl;
      G4cout << "Incoming photon energy: " << photonEnergy0/keV << " keV" << G4endl;
      G4cout << "-----------------------------------------------------------" << G4endl;
      G4cout << "Scattered photon: " << photonEnergy1/keV << " keV" << G4endl;
      if (electron)
	electronEnergy = eKineticEnergy;
      G4cout << "Scattered electron " << electronEnergy/keV << " keV" << G4endl;
      G4cout << "Fluorescence: " << energyInFluorescence/keV << " keV" << G4endl;
      G4cout << "Local energy deposit " << localEnergyDeposit/keV << " keV" << G4endl;
      G4cout << "Total final state: " << (photonEnergy1+electronEnergy+energyInFluorescence+
					  localEnergyDeposit)/keV << 
	" keV" << G4endl;
      G4cout << "-----------------------------------------------------------" << G4endl;
    }
  if (verboseLevel > 0)
    {
      G4double energyDiff = std::fabs(photonEnergy1+
				      electronEnergy+energyInFluorescence+
				      localEnergyDeposit-photonEnergy0);
      if (energyDiff > 0.05*keV)
	G4cout << "Warning from G4Penelope08Compton: problem with energy conservation: " << 
	  (photonEnergy1+electronEnergy+energyInFluorescence+localEnergyDeposit)/keV << 
	  " keV (final) vs. " << 
	  photonEnergy0/keV << " keV (initial)" << G4endl;
    }    
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....

G4double G4Penelope08ComptonModel::DifferentialCrossSection(G4double cosTheta,G4double energy,
							    G4PenelopeOscillator* osc)
{
  //
  // Penelope model. 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(); 

  const G4double k2 = std::sqrt(2.);
  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*std::exp(0.5-(k1+k2*x)*(k1+k2*x));    
  else    
    sia = 0.5*std::exp(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....

void G4Penelope08ComptonModel::ActivateAuger(G4bool augerbool)
{
  if (!DeexcitationFlag() && augerbool)
    {
      G4cout << "WARNING - G4Penelope08ComptonModel" << G4endl;
      G4cout << "The use of the Atomic Deexcitation Manager is set to false " << G4endl;
      G4cout << "Therefore, Auger electrons will be not generated anyway" << G4endl;
    }
  deexcitationManager.ActivateAugerElectronProduction(augerbool);
  if (verboseLevel > 1)
    G4cout << "Auger production set to " << augerbool << G4endl;
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....

G4double G4Penelope08ComptonModel::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 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 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 d = Weights[0]*
	(DifferentialCrossSection(b+c,energy,osc)+DifferentialCrossSection(b-c,energy,osc));
      
      for (G4int j=1;j<npoints;j++)
	{
	  c=a*Abscissas[j];
	  d += Weights[j]*
	    (DifferentialCrossSection(b+c,energy,osc)+DifferentialCrossSection(b-c,energy,osc));
	}    
      G4double s1=d*a;
      a=0.5*(xc-xb);
      b=0.5*(xc+xb);
      c=a*Abscissas[0];
      d=Weights[0]*
	(DifferentialCrossSection(b+c,energy,osc)+DifferentialCrossSection(b-c,energy,osc));
      
      for (G4int j=1;j<npoints;j++)
	{
	  c=a*Abscissas[j];
	  d += Weights[j]*
	    (DifferentialCrossSection(b+c,energy,osc)+DifferentialCrossSection(b-c,energy,osc));
	}    
      G4double s2=d*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 << "G4Penelope08ComptonModel: " << 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 G4Penelope08ComptonModel::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*std::log(t0)-(1.0/t0);

  G4PenelopeOscillatorTable* theTable = oscManager->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*std::log(tau)-(1.0/tau);
	  G4double stre = theOsc->GetOscillatorStrength();

	  cs += stre*(csu-csl);
	}
    }

  cs /= (ek*eks);

  return cs;

}