//
// ********************************************************************
// * 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: G4Penelope08RayleighModel.cc,v 1.3 2010/07/28 07:09:16 pandola Exp $
// GEANT4 tag $Name: geant4-09-04-ref-00 $
//
// Author: Luciano Pandola
//
// History:
// --------
// 03 Dec 2009   L Pandola    First implementation
//

#include "G4Penelope08RayleighModel.hh"
#include "G4PenelopeSamplingData.hh"
#include "G4ParticleDefinition.hh"
#include "G4MaterialCutsCouple.hh"
#include "G4ProductionCutsTable.hh"
#include "G4DynamicParticle.hh"
#include "G4PhysicsTable.hh"
#include "G4ElementTable.hh"
#include "G4Element.hh"
#include "G4PhysicsFreeVector.hh"


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


G4Penelope08RayleighModel::G4Penelope08RayleighModel(const G4ParticleDefinition*,
                                             const G4String& nam)
  :G4VEmModel(nam),isInitialised(false),logAtomicCrossSection(0),   
   atomicFormFactor(0),logFormFactorTable(0),pMaxTable(0),samplingTable(0)
{
  fIntrinsicLowEnergyLimit = 100.0*eV;
  fIntrinsicHighEnergyLimit = 100.0*GeV;
  //  SetLowEnergyLimit(fIntrinsicLowEnergyLimit);
  SetHighEnergyLimit(fIntrinsicHighEnergyLimit);
  //
  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 

  //build the energy grid. It is the same for all materials
  G4double logenergy = std::log(fIntrinsicLowEnergyLimit/2.);
  G4double logmaxenergy = std::log(1.5*fIntrinsicHighEnergyLimit);
  //finer grid below 160 keV
  G4double logtransitionenergy = std::log(160*keV); 
  G4double logfactor1 = std::log(10.)/250.;
  G4double logfactor2 = logfactor1*10;
  logEnergyGridPMax.push_back(logenergy);
  do{
    if (logenergy < logtransitionenergy)
      logenergy += logfactor1;
    else
      logenergy += logfactor2;
    logEnergyGridPMax.push_back(logenergy);      
  }while (logenergy < logmaxenergy);
}

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

G4Penelope08RayleighModel::~G4Penelope08RayleighModel()
{
  std::map <const G4int,G4PhysicsFreeVector*>::iterator i;
  if (logAtomicCrossSection)
    {
      for (i=logAtomicCrossSection->begin();i != logAtomicCrossSection->end();i++)
	if (i->second) delete i->second;
      delete logAtomicCrossSection;
     }

   if (atomicFormFactor)
     {
       for (i=atomicFormFactor->begin();i != atomicFormFactor->end();i++)
	 if (i->second) delete i->second;
       delete atomicFormFactor;
     }

  ClearTables();
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
void G4Penelope08RayleighModel::ClearTables()
{
   std::map <const G4Material*,G4PhysicsFreeVector*>::iterator i;
 
   if (logFormFactorTable)
     {
       for (i=logFormFactorTable->begin(); i != logFormFactorTable->end(); i++)
	 if (i->second) delete i->second;
       delete logFormFactorTable;
       logFormFactorTable = 0; //zero explicitely
     }

   if (pMaxTable)
     {
       for (i=pMaxTable->begin(); i != pMaxTable->end(); i++)
	 if (i->second) delete i->second;
       delete pMaxTable;
       pMaxTable = 0; //zero explicitely
     }

   std::map<const G4Material*,G4PenelopeSamplingData*>::iterator ii;
   if (samplingTable)
     {
       for (ii=samplingTable->begin(); ii != samplingTable->end(); ii++)
	 if (ii->second) delete ii->second;
       delete samplingTable;
       samplingTable = 0; //zero explicitely
     }     

   return;
}

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

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

  //clear tables depending on materials, not the atomic ones
  ClearTables();
  
  //create new tables
  //
  // logAtomicCrossSection and atomicFormFactor are created only once,
  // since they are never cleared
  if (!logAtomicCrossSection)
    logAtomicCrossSection = new std::map<const G4int,G4PhysicsFreeVector*>;
  if (!atomicFormFactor)
    atomicFormFactor = new std::map<const G4int,G4PhysicsFreeVector*>;

  if (!logFormFactorTable)
    logFormFactorTable = new std::map<const G4Material*,G4PhysicsFreeVector*>;
  if (!pMaxTable)
    pMaxTable = new std::map<const G4Material*,G4PhysicsFreeVector*>;
  if (!samplingTable)
    samplingTable = new std::map<const G4Material*,G4PenelopeSamplingData*>;


  if (verboseLevel > 0) {
    G4cout << "Penelope08 Rayleigh 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 G4Penelope08RayleighModel::ComputeCrossSectionPerAtom(const G4ParticleDefinition*,
							       G4double energy,
							       G4double Z,
							       G4double,
							       G4double,
							       G4double)
{
  // Cross section of Rayleigh scattering in Penelope2008 is calculated by the EPDL97 
  // tabulation, Cuellen et al. (1997), with non-relativistic form factors from Hubbel 
  // et al. J. Phys. Chem. Ref. Data 4 (1975) 471; Erratum ibid. 6 (1977) 615.

   if (verboseLevel > 3)
    G4cout << "Calling CrossSectionPerAtom() of G4Penelope08RayleighModel" << G4endl;
 
   G4int iZ = (G4int) Z;

   //read data files
   if (!logAtomicCrossSection->count(iZ))
     ReadDataFile(iZ);
   //now it should be ok
   if (!logAtomicCrossSection->count(iZ))
     {
       G4cout << "Problem in G4Penelope08RayleighModel::ComputeCrossSectionPerAtom" 
	      << G4endl;
       G4Exception();
     }

   G4double cross = 0;

   G4PhysicsFreeVector* atom = logAtomicCrossSection->find(iZ)->second;
   if (!atom)
     {
       G4cout << "Problem in G4Penelope08RayleighModel::ComputeCrossSectionPerAtom" 
	 << G4endl;
       G4Exception();
     }
   G4double logene = std::log(energy);
   G4double logXS = atom->Value(logene);
   cross = std::exp(logXS);

   if (verboseLevel > 2)
    G4cout << "Rayleigh cross section at " << energy/keV << " keV for Z=" << Z << 
      " = " << cross/barn << " barn" << G4endl;
    return cross;
}


//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
void G4Penelope08RayleighModel::BuildFormFactorTable(const G4Material* material)
{
 
  /*
    1) get composition and equivalent molecular density
  */
  
  G4int nElements = material->GetNumberOfElements();
  const G4ElementVector* elementVector = material->GetElementVector();
  const G4double* fractionVector = material->GetFractionVector();

  std::vector<G4double> *StechiometricFactors = new std::vector<G4double>;
  for (G4int i=0;i<nElements;i++)
    {
      G4double fraction = fractionVector[i];
      G4double atomicWeigth = (*elementVector)[i]->GetA()/(g/mole);
      StechiometricFactors->push_back(fraction/atomicWeigth);
    }
  //Find max
  G4double MaxStechiometricFactor = 0.;
  for (G4int i=0;i<nElements;i++)
    {
      if ((*StechiometricFactors)[i] > MaxStechiometricFactor)
        MaxStechiometricFactor = (*StechiometricFactors)[i];
    }
  if (MaxStechiometricFactor<1e-16)
    {
      G4cout << "G4Penelope08RayleighModel::BuildFormFactorTable" << G4endl;
      G4cout << "Problem with the mass composition of " << material->GetName() << G4endl;
      G4Exception();
    }
  //Normalize
  for (G4int i=0;i<nElements;i++)
    (*StechiometricFactors)[i] /=  MaxStechiometricFactor;
 
  // Equivalent atoms per molecule
  G4double atomsPerMolecule = 0;
  for (G4int i=0;i<nElements;i++)
    atomsPerMolecule += (*StechiometricFactors)[i]; 
 
  /*
    CREATE THE FORM FACTOR TABLE
  */
  G4PhysicsFreeVector* theFFVec = new G4PhysicsFreeVector(logQSquareGrid.size());
  theFFVec->SetSpline(true);

  for (size_t k=0;k<logQSquareGrid.size();k++)
    {
      G4double ff2 = 0; //squared form factor
      for (G4int i=0;i<nElements;i++)
	{
	  G4int iZ = (G4int) (*elementVector)[i]->GetZ();
	  G4PhysicsFreeVector* theAtomVec = atomicFormFactor->find(iZ)->second;
	  G4double f = (*theAtomVec)[k]; //the q-grid is always the same	    
	  ff2 += f*f*(*StechiometricFactors)[i];
	}
      if (ff2)
	theFFVec->PutValue(k,logQSquareGrid[k],std::log(ff2)); //NOTICE: THIS IS log(Q^2) vs. log(F^2)
    }
  logFormFactorTable->insert(std::make_pair(material,theFFVec));

  delete StechiometricFactors;
  
  return;
}


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

void G4Penelope08RayleighModel::SampleSecondaries(std::vector<G4DynamicParticle*>* ,
						const G4MaterialCutsCouple* couple,
						const G4DynamicParticle* aDynamicGamma,
						G4double,
						G4double)
{
  // Sampling of the Rayleigh final state (namely, scattering angle of the photon) 
  // from the Penelope2008 model. The scattering angle is sampled from the atomic 
  // cross section dOmega/d(cosTheta) from Born ("Atomic Phyisics", 1969), disregarding 
  // anomalous scattering effects. The Form Factor F(Q) function which appears in the 
  // analytical cross section is retrieved via the method GetFSquared(); atomic data 
  // are tabulated for F(Q). Form factor for compounds is calculated according to 
  // the additivity rule. The sampling from the F(Q) is made via a Rational Inverse 
  // Transform with Aliasing (RITA) algorithm; RITA parameters are calculated once 
  // for each material and managed by G4PenelopeSamplingData objects.
  // The sampling algorithm (rejection method) has efficiency 67% at low energy, and 
  // increases with energy. For E=100 keV the efficiency is 100% and 86% for 
  // hydrogen and uranium, respectively.

  if (verboseLevel > 3)
    G4cout << "Calling SamplingSecondaries() of G4Penelope08RayleighModel" << G4endl;

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

  G4ParticleMomentum photonDirection0 = aDynamicGamma->GetMomentumDirection();
  
  const G4Material* theMat = couple->GetMaterial();
  
  //1) Verify if tables are ready
  if (!pMaxTable || !samplingTable)
    {
      G4cout << " G4Penelope08RayleighModel::SampleSecondaries" << G4endl;
      G4cout << "Looks like the model is _not_ properly initialized" << G4endl;
      G4Exception();
    }
  
  //2) retrieve or build the sampling table
  if (!(samplingTable->count(theMat)))
    InitializeSamplingAlgorithm(theMat);
  G4PenelopeSamplingData* theDataTable = samplingTable->find(theMat)->second;
  
  //3) retrieve or build the pMax data
  if (!pMaxTable->count(theMat))
    GetPMaxTable(theMat);
  G4PhysicsFreeVector* thePMax = pMaxTable->find(theMat)->second;

  G4double cosTheta = 1.0;
  
  //OK, ready to go!
  G4double qmax = 2.0*photonEnergy0/electron_mass_c2; //this is non-dimensional now

  if (qmax < 1e-10) //very low momentum transfer
    {
      G4bool loopAgain=false;
      do
	{
	  loopAgain = false;
	  cosTheta = 1.0-2.0*G4UniformRand();
	  G4double G = 0.5*(1+cosTheta*cosTheta);
	  if (G4UniformRand()>G)
	    loopAgain = true;
	}while(loopAgain);
    }
  else //larger momentum transfer
    {
      size_t nData = theDataTable->GetNumberOfStoredPoints();
      G4double LastQ2inTheTable = theDataTable->GetX(nData-1);
      G4double q2max = std::min(qmax*qmax,LastQ2inTheTable);

      G4bool loopAgain = false;
      G4double MaxPValue = thePMax->Value(photonEnergy0);
      G4double xx=0;
      
      //Sampling by rejection method. The rejection function is 
      //G = 0.5*(1+cos^2(theta))
      //
      do{
	loopAgain = false;
	G4double RandomMax = G4UniformRand()*MaxPValue;
	xx = theDataTable->SampleValue(RandomMax);
	//xx is a random value of q^2 in (0,q2max),sampled according to 
	//F(Q^2) via the RITA algorithm
	if (xx > q2max)
	  loopAgain = true;
	cosTheta = 1.0-2.0*xx/q2max;
	G4double G = 0.5*(1+cosTheta*cosTheta);
	if (G4UniformRand()>G)
	  loopAgain = true;
      }while(loopAgain);
    }
  
  G4double sinTheta = std::sqrt(1-cosTheta*cosTheta);
 
  // Scattered photon angles. ( Z - axis along the parent photon)
  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) ;
  fParticleChange->SetProposedKineticEnergy(photonEnergy0) ;
 
  return;
}


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

void G4Penelope08RayleighModel::ReadDataFile(const G4int Z)
{
  if (verboseLevel > 2)
    {
      G4cout << "G4Penelope08RayleighModel::ReadDataFile()" << G4endl;
      G4cout << "Going to read Rayleigh data files for Z=" << Z << G4endl;
    }

  char* path = getenv("G4LEDATA");
  if (!path)
    {
      G4String excep = "G4Penelope08RayleighModel - G4LEDATA environment variable not set!";
      G4Exception(excep);
    }

  /*
    Read first the cross section file
  */
  std::ostringstream ost;
  if (Z>9)
    ost << path << "/penelope/rayleigh/pdgra" << Z << ".p08";
  else
    ost << path << "/penelope/rayleigh/pdgra0" << Z << ".p08";
  std::ifstream file(ost.str().c_str());
  if (!file.is_open())
    {
      G4String excep = "G4Penelope08RayleighModel - data file " + G4String(ost.str()) + " not found!";
      G4Exception(excep);
    }
  G4int readZ =0;
  size_t nPoints= 0;
  file >> readZ >> nPoints;
  //check the right file is opened.
  if (readZ != Z || nPoints <= 0)
    {
      G4cout << "G4Penelope08RayleighModel::ReadDataFile()" << G4endl;
      G4cout << "Corrupted data file for Z=" << Z << G4endl;
      G4Exception();
    }  
  G4PhysicsFreeVector* theVec = new G4PhysicsFreeVector((size_t)nPoints);
  G4double ene=0,f1=0,f2=0,xs=0;
  for (size_t i=0;i<nPoints;i++)
    {
      file >> ene >> f1 >> f2 >> xs;
      //dimensional quantities
      ene *= eV;
      xs *= cm2;
      theVec->PutValue(i,std::log(ene),std::log(xs));
      if (file.eof() && i != (nPoints-1)) //file ended too early
	{
	  G4cout << "G4Penelope08RayleighModel::ReadDataFile()" << G4endl;
	  G4cout << "Corrupted data file for Z=" << Z << G4endl;
	  G4cout << "Found less than " << nPoints << "entries " <<G4endl;
	  G4Exception();
	}
    }
  if (!logAtomicCrossSection)
    {
      G4cout << "G4Penelope08RayleighModel::ReadDataFile()" << G4endl;
      G4cout << "Problem with allocation of logAtomicCrossSection data table " << G4endl;
      G4Exception();
    }
  logAtomicCrossSection->insert(std::make_pair(Z,theVec));
  file.close();

  /*
    Then read the form factor file
  */
  std::ostringstream ost2;
  if (Z>9)
    ost2 << path << "/penelope/rayleigh/pdaff" << Z << ".p08";
  else
    ost2 << path << "/penelope/rayleigh/pdaff0" << Z << ".p08";
  file.open(ost2.str().c_str());
  if (!file.is_open())
    {
      G4String excep = "G4Penelope08RayleighModel - data file " + G4String(ost2.str()) + " not found!";
      G4Exception(excep);
    }
  file >> readZ >> nPoints;
  //check the right file is opened.
  if (readZ != Z || nPoints <= 0)
    {
      G4cout << "G4Penelope08RayleighModel::ReadDataFile()" << G4endl;
      G4cout << "Corrupted data file for Z=" << Z << G4endl;
      G4Exception();
    }  
  G4PhysicsFreeVector* theFFVec = new G4PhysicsFreeVector((size_t)nPoints);
  G4double q=0,ff=0,incoh=0;
  G4bool fillQGrid = false;
  //fill this vector only the first time.
  if (!logQSquareGrid.size())
    fillQGrid = true;
  for (size_t i=0;i<nPoints;i++)
    {
      file >> q >> ff >> incoh;
      //q and ff are dimensionless (q is in units of (m_e*c)
      theFFVec->PutValue(i,q,ff);
      if (fillQGrid)
	{
	  logQSquareGrid.push_back(2.0*std::log(q));
	}
      if (file.eof() && i != (nPoints-1)) //file ended too early
	{
	  G4cout << "G4Penelope08RayleighModel::ReadDataFile()" << G4endl;
	  G4cout << "Corrupted data file for Z=" << Z << G4endl;
	  G4cout << "Found less than " << nPoints << "entries " <<G4endl;
	  G4Exception();
	}
    }
  if (!atomicFormFactor)
    {
      G4cout << "G4Penelope08RayleighModel::ReadDataFile()" << G4endl;
      G4cout << "Problem with allocation of atomicFormFactor data table " << G4endl;
      G4Exception();
    }
  atomicFormFactor->insert(std::make_pair(Z,theFFVec));
  file.close();
  return;
}

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

G4double G4Penelope08RayleighModel::GetFSquared(const G4Material* mat, const G4double QSquared)
{
  G4double f2 = 0;
  G4double logQSquared = std::log(QSquared);
  //last value of the table
  G4double maxlogQ2 = logQSquareGrid[logQSquareGrid.size()-1];
  //If the table has not been built for the material, do it!
  if (!logFormFactorTable->count(mat))
    BuildFormFactorTable(mat);
  
  //now it should  be all right
  G4PhysicsFreeVector* theVec = logFormFactorTable->find(mat)->second;

  if (!theVec)
    {
      G4cout << " G4Penelope08RayleighModel::GetFSquared()" << G4endl;
      G4cout << "Unable to retrieve F squared table for " << mat->GetName();
      G4Exception();
    }
  if (logQSquared < -20) // Q < 1e-9
    {
      G4double logf2 = (*theVec)[0]; //first value of the table
      f2 = std::exp(logf2);
    }
  else if (logQSquared > maxlogQ2)
    f2 =0;
  else
    {
      //log(Q^2) vs. log(F^2)
      G4double logf2 = theVec->Value(logQSquared);
      f2 = std::exp(logf2);

    }
  if (verboseLevel > 3)
    {
      G4cout << "G4Penelope08RayleighModel::GetFSquared() in " << mat->GetName() << G4endl;  
      G4cout << "Q^2 = " <<  QSquared << " (units of 1/(m_e*c); F^2 = " << f2 << G4endl;
    }
  return f2;
}

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

void G4Penelope08RayleighModel::InitializeSamplingAlgorithm(const G4Material* mat)
{
  G4double q2min = 0;
  G4double q2max = 0;
  const size_t np = 150; //hard-coded in Penelope
  for (size_t i=1;i<logQSquareGrid.size();i++)
    {
      G4double Q2 = std::exp(logQSquareGrid[i]);
      if (GetFSquared(mat,Q2) >  1e-35)
	{
	  q2max = std::exp(logQSquareGrid[i-1]);
	}
    }
  
  size_t nReducedPoints = np/4;

  //check for errors
  if (np < 16)
    {
      G4cout << "G4Penelope08RayleighModel::InitializeSamplingAlgorithm" << G4endl;
      G4cout << "Too few points to initialize the sampling algorithm" << G4endl;
      G4Exception();
    }
  if (q2min > (q2max-1e-10))
    {
      G4cout << "G4Penelope08RayleighModel::InitializeSamplingAlgorithm" << G4endl;
      G4cout << "Too narrow grid to initialize the sampling algorithm" << G4endl;
      G4Exception();
    }

  //This is subroutine RITAI0 of Penelope
  //Create an object of type G4Penelope08RayleighSamplingData --> store in a map::Material*

  //temporary vectors --> Then everything is stored in G4PenelopeSamplingData
  G4DataVector* x = new G4DataVector();
  
  /*******************************************************************************
    Start with a grid of NUNIF points uniformly spaced in the interval q2min,q2max
  ********************************************************************************/
  size_t NUNIF = std::min(std::max(((size_t)8),nReducedPoints),np/2);
  const G4int nip = 51; //hard-coded in Penelope

  G4double dx = (q2max-q2min)/((G4double) NUNIF-1);
  x->push_back(q2min);
  for (size_t i=1;i<NUNIF-1;i++)
    {
      G4double app = q2min + i*dx;
      x->push_back(app); //increase 
    }
  x->push_back(q2max);
  
  if (verboseLevel> 3)
    G4cout << "Vector x has " << x->size() << " points, while NUNIF = " << NUNIF << G4endl;

  //Allocate temporary storage vectors
  G4DataVector* area = new G4DataVector();
  G4DataVector* a = new G4DataVector();
  G4DataVector* b = new G4DataVector();
  G4DataVector* c = new G4DataVector();
  G4DataVector* err = new G4DataVector();

  for (size_t i=0;i<NUNIF-1;i++) //build all points but the last
    {      
      //Temporary vectors for this loop
      G4DataVector* pdfi = new G4DataVector();
      G4DataVector* pdfih = new G4DataVector();
      G4DataVector* sumi = new G4DataVector();

      G4double dxi = ((*x)[i+1]-(*x)[i])/(G4double (nip-1));
      G4double pdfmax = 0;
      for (G4int k=0;k<nip;k++)
	{
	  G4double xik = (*x)[i]+k*dxi;
	  G4double pdfk = std::max(GetFSquared(mat,xik),0.);
	  pdfi->push_back(pdfk);
	  pdfmax = std::max(pdfmax,pdfk);	
	  if (k < (nip-1))
	    {
	      G4double xih = xik + 0.5*dxi;
	      G4double pdfIK = std::max(GetFSquared(mat,xih),0.);
	      pdfih->push_back(pdfIK);
	      pdfmax = std::max(pdfmax,pdfIK);
	    }
	}
    
      //Simpson's integration
      G4double cons = dxi*0.5*(1./3.);
      sumi->push_back(0.);
      for (G4int k=1;k<nip;k++)
	{
	  G4double previous = (*sumi)[k-1];
	  G4double next = previous + cons*((*pdfi)[k-1]+4.0*(*pdfih)[k-1]+(*pdfi)[k]);
	  sumi->push_back(next);
	}
    
      G4double lastIntegral = (*sumi)[sumi->size()-1];
      area->push_back(lastIntegral);
      //Normalize cumulative function
      G4double factor = 1.0/lastIntegral;
      for (size_t k=0;k<sumi->size();k++)
	(*sumi)[k] *= factor;
      
      //When the PDF vanishes at one of the interval end points, its value is modified
      if ((*pdfi)[0] < 1e-35) 
	(*pdfi)[0] = 1e-5*pdfmax;
      if ((*pdfi)[pdfi->size()-1] < 1e-35)
	(*pdfi)[pdfi->size()-1] = 1e-5*pdfmax;

      G4double pli = (*pdfi)[0]*factor;
      G4double pui = (*pdfi)[pdfi->size()-1]*factor;
      G4double B_temp = 1.0-1.0/(pli*pui*dx*dx);
      G4double A_temp = (1.0/(pli*dx))-1.0-B_temp;
      G4double C_temp = 1.0+A_temp+B_temp;
      if (C_temp < 1e-35)
	{
	  a->push_back(0.);
	  b->push_back(0.);
	  c->push_back(1.);	  
	}
      else
	{
	  a->push_back(A_temp);
	  b->push_back(B_temp);
	  c->push_back(C_temp);
	}

      //OK, now get ERR(I), the integral of the absolute difference between the rational interpolation 
      //and the true pdf, extended over the interval (X(I),X(I+1))
      G4int icase = 1; //loop code
      G4bool reLoop = false;
      err->push_back(0.);
      do
	{
	  reLoop = false;
	  (*err)[i] = 0.; //zero variable
	  for (G4int k=0;k<nip;k++)
	    {
	      G4double rr = (*sumi)[k];
	      G4double pap = (*area)[i]*(1.0+((*a)[i]+(*b)[i]*rr)*rr)*(1.0+((*a)[i]+(*b)[i]*rr)*rr)/
		((1.0-(*b)[i]*rr*rr)*(*c)[i]*((*x)[i+1]-(*x)[i]));
	      if (k == 0 || k == nip-1)
		(*err)[i] += 0.5*std::fabs(pap-(*pdfi)[k]);
	      else
		(*err)[i] += std::fabs(pap-(*pdfi)[k]);
	    }
	  (*err)[i] *= dxi;
      
	  //If err(I) is too large, the pdf is approximated by a uniform distribution
	  if ((*err)[i] > 0.1*(*area)[i] && icase == 1) 
	    {
	      (*b)[i] = 0;
	      (*a)[i] = 0;
	      (*c)[i] = 1.;
	      icase = 2;
	      reLoop = true;
	    }
	}while(reLoop);

      if (pdfi) delete pdfi;
      if (pdfih) delete pdfih;
      if (sumi) delete sumi;
    } //end of first loop over i

  //Now assign last point
  (*x)[x->size()-1] = q2max;
  a->push_back(0.);
  b->push_back(0.);
  c->push_back(0.);
  err->push_back(0.);
  area->push_back(0.);

  if (x->size() != NUNIF || a->size() != NUNIF || 
      err->size() != NUNIF || area->size() != NUNIF)
    {
      G4cout << "G4Penelope08RayleighModel::InitializeSamplingAlgorithm" << G4endl;
      G4cout << "Problem in building the Table for Sampling: array dimensions do not match " << G4endl;
      G4Exception();
    }
  
  /*******************************************************************************
   New grid points are added by halving the sub-intervals with the largest absolute error
  This is done up to np=150 points in the grid
  ********************************************************************************/
  do
    {
      G4double maxError = 0.0;
      size_t iErrMax = 0;
      for (size_t i=0;i<err->size()-2;i++) 
	{
	  //maxError is the lagest of the interval errors err[i]
	  if ((*err)[i] > maxError)
	    {
	      maxError = (*err)[i];
	      iErrMax = i;
	    }
	}
      
      //OK, now I have to insert one new point in the position iErrMax
      G4double newx = 0.5*((*x)[iErrMax]+(*x)[iErrMax+1]);
      
      x->insert(x->begin()+iErrMax+1,newx);
      //Add place-holders in the other vectors
      area->insert(area->begin()+iErrMax+1,0.);
      a->insert(a->begin()+iErrMax+1,0.);
      b->insert(b->begin()+iErrMax+1,0.);
      c->insert(c->begin()+iErrMax+1,0.);
      err->insert(err->begin()+iErrMax+1,0.);
        
      //Now calculate the other parameters
      for (size_t i=iErrMax;i<=iErrMax+1;i++)
	{
	  //Temporary vectors for this loop
	  G4DataVector* pdfi = new G4DataVector();
	  G4DataVector* pdfih = new G4DataVector();
	  G4DataVector* sumi = new G4DataVector();
	  
	  G4double dx = (*x)[i+1]-(*x)[i];
	  G4double dxi = ((*x)[i+1]-(*x)[i])/(G4double (nip-1));
	  G4double pdfmax = 0;
	  for (G4int k=0;k<nip;k++)
	    {
	      G4double xik = (*x)[i]+k*dxi;
	      G4double pdfk = std::max(GetFSquared(mat,xik),0.);
	      pdfi->push_back(pdfk);
	      pdfmax = std::max(pdfmax,pdfk);	
	      if (k < (nip-1))
		{
		  G4double xih = xik + 0.5*dxi;
		  G4double pdfIK = std::max(GetFSquared(mat,xih),0.);
		  pdfih->push_back(pdfIK);
		  pdfmax = std::max(pdfmax,pdfIK);
		}
	    }
	  
	  //Simpson's integration
	  G4double cons = dxi*0.5*(1./3.);
	  sumi->push_back(0.);
	  for (G4int k=1;k<nip;k++)
	    {
	      G4double previous = (*sumi)[k-1];
	      G4double next = previous + cons*((*pdfi)[k-1]+4.0*(*pdfih)[k-1]+(*pdfi)[k]);
	      sumi->push_back(next);
	    }
	  G4double lastIntegral = (*sumi)[sumi->size()-1];
	  (*area)[i] = lastIntegral;
	  
	  //Normalize cumulative function
	  G4double factor = 1.0/lastIntegral;
	  for (size_t k=0;k<sumi->size();k++)
	    (*sumi)[k] *= factor;
	  
	  //When the PDF vanishes at one of the interval end points, its value is modified
	  if ((*pdfi)[0] < 1e-35) 
	    (*pdfi)[0] = 1e-5*pdfmax;
	  if ((*pdfi)[pdfi->size()-1] < 1e-35)
	    (*pdfi)[pdfi->size()-1] = 1e-5*pdfmax;
	  
	  G4double pli = (*pdfi)[0]*factor;
	  G4double pui = (*pdfi)[pdfi->size()-1]*factor;
	  G4double B_temp = 1.0-1.0/(pli*pui*dx*dx);
	  G4double A_temp = (1.0/(pli*dx))-1.0-B_temp;
	  G4double C_temp = 1.0+A_temp+B_temp;
	  if (C_temp < 1e-35)
	    {
	      (*a)[i]= 0.;
	      (*b)[i] = 0.;
	      (*c)[i] = 0;
	    }
	  else
	    {
	      (*a)[i]= A_temp;
	      (*b)[i] = B_temp;
	      (*c)[i] = C_temp;
	    }
	  //OK, now get ERR(I), the integral of the absolute difference between the rational interpolation 
	  //and the true pdf, extended over the interval (X(I),X(I+1))
	  G4int icase = 1; //loop code
	  G4bool reLoop = false;
	  do
	    {
	      reLoop = false;
	      (*err)[i] = 0.; //zero variable
	      for (G4int k=0;k<nip;k++)
		{
		  G4double rr = (*sumi)[k];
		  G4double pap = (*area)[i]*(1.0+((*a)[i]+(*b)[i]*rr)*rr)*(1.0+((*a)[i]+(*b)[i]*rr)*rr)/
		    ((1.0-(*b)[i]*rr*rr)*(*c)[i]*((*x)[i+1]-(*x)[i]));
		  if (k == 0 || k == nip-1)
		    (*err)[i] += 0.5*std::fabs(pap-(*pdfi)[k]);
		  else
		    (*err)[i] += std::fabs(pap-(*pdfi)[k]);
		}
	      (*err)[i] *= dxi;
	      
	      //If err(I) is too large, the pdf is approximated by a uniform distribution
	      if ((*err)[i] > 0.1*(*area)[i] && icase == 1) 
		{
		  (*b)[i] = 0;
		  (*a)[i] = 0;
		  (*c)[i] = 1.;
		  icase = 2;
		  reLoop = true;
		}
	    }while(reLoop);
	  if (pdfi) delete pdfi;
	  if (pdfih) delete pdfih;
	  if (sumi) delete sumi;
	}
    }while(x->size() < np);

  if (x->size() != np || a->size() != np || 
      err->size() != np || area->size() != np)
    {
      G4cout << "G4Penelope08RayleighModel::InitializeSamplingAlgorithm" << G4endl;
      G4cout << "Problem in building the extended Table for Sampling: array dimensions do not match " << G4endl;
      G4Exception();
    }

  /*******************************************************************************
   Renormalization
  ********************************************************************************/
  G4double ws = 0;
  for (size_t i=0;i<np-1;i++)
    ws += (*area)[i];
  ws = 1.0/ws;
  G4double errMax = 0;
  for (size_t i=0;i<np-1;i++)
    {
      (*area)[i] *= ws;
      (*err)[i] *= ws;
      errMax = std::max(errMax,(*err)[i]);
    }

  //Vector with the normalized cumulative distribution
  G4DataVector* PAC = new G4DataVector();
  PAC->push_back(0.);
  for (size_t i=0;i<np-1;i++)
    {
      G4double previous = (*PAC)[i];
      PAC->push_back(previous+(*area)[i]);
    }
  (*PAC)[PAC->size()-1] = 1.;
		     
  /*******************************************************************************
  Pre-calculated limits for the initial binary search for subsequent sampling
  ********************************************************************************/

  //G4DataVector* ITTL = new G4DataVector();
  std::vector<size_t> *ITTL = new std::vector<size_t>;
  std::vector<size_t> *ITTU = new std::vector<size_t>;

  //Just create place-holders
  for (size_t i=0;i<np;i++)
    {
      ITTL->push_back(0);
      ITTU->push_back(0);
    }

  G4double bin = 1.0/(np-1);
  (*ITTL)[0]=0;
  for (size_t i=1;i<(np-1);i++)
    {
      G4double ptst = i*bin; 
      G4bool found = false;
      for (size_t j=(*ITTL)[i-1];j<np && !found;j++)
	{
	  if ((*PAC)[j] > ptst)
	    {
	      (*ITTL)[i] = j-1;
	      (*ITTU)[i-1] = j;
	      found = true;
	    }
	}
    }
  (*ITTU)[ITTU->size()-2] = ITTU->size()-1;
  (*ITTU)[ITTU->size()-1] = ITTU->size()-1;
  (*ITTL)[ITTL->size()-1] = ITTU->size()-2;

  if (ITTU->size() != np || ITTU->size() != np)
    {
      G4cout << "G4Penelope08RayleighModel::InitializeSamplingAlgorithm" << G4endl;
      G4cout << "Problem in building the Limit Tables for Sampling: array dimensions do not match " << G4endl;
      G4Exception();
    }


  /********************************************************************************
    Copy tables
  ********************************************************************************/
  G4PenelopeSamplingData* theTable = new G4PenelopeSamplingData(np);
  for (size_t i=0;i<np;i++)
    {
      theTable->AddPoint((*x)[i],(*PAC)[i],(*a)[i],(*b)[i],(*ITTL)[i],(*ITTU)[i]);
    }

  if (verboseLevel > 2)
    {
      G4cout << "*************************************************************************" << G4endl;
      G4cout << "Sampling table for Penelope Rayleigh scattering in " << mat->GetName() << G4endl;
      theTable->DumpTable();
    } 
  samplingTable->insert(std::make_pair(mat,theTable));

 
  //Clean up temporary vectors
  if (x) delete x;
  if (a) delete a;
  if (b) delete b;
  if (c) delete c;
  if (err) delete err;
  if (area) delete area;
  if (PAC) delete PAC;
  if (ITTL) delete ITTL;
  if (ITTU) delete ITTU;

  //DONE!
  return;
  
}

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

void G4Penelope08RayleighModel::GetPMaxTable(const G4Material* mat)
{
  if (!pMaxTable)
    {
      G4cout << "G4Penelope08RayleighModel::BuildPMaxTable" << G4endl;
      G4cout << "Going to instanziate the pMaxTable !" << G4endl;
      G4cout << "That should _not_ be here! " << G4endl;
      pMaxTable = new std::map<const G4Material*,G4PhysicsFreeVector*>;
    }
  //check if the table is already there
  if (pMaxTable->count(mat))
    return;

  //otherwise build it
  if (!samplingTable)
    {
      G4cout << "G4Penelope08RayleighModel::BuildPMaxTable" << G4endl;
      G4cout << "WARNING: samplingTable is not properly instantiated" << G4endl;
      G4Exception();
    }

  if (!samplingTable->count(mat))
    InitializeSamplingAlgorithm(mat);      
  
  G4PenelopeSamplingData *theTable = samplingTable->find(mat)->second;
  size_t tablePoints = theTable->GetNumberOfStoredPoints();

  size_t nOfEnergyPoints = logEnergyGridPMax.size();
  G4PhysicsFreeVector* theVec = new G4PhysicsFreeVector(nOfEnergyPoints);

  const size_t nip = 51; //hard-coded in Penelope

  for (size_t ie=0;ie<logEnergyGridPMax.size();ie++)
    {
      G4double energy = std::exp(logEnergyGridPMax[ie]);
      G4double Qm = 2.0*energy/electron_mass_c2; //this is non-dimensional now
      G4double Qm2 = Qm*Qm;
      G4double firstQ2 = theTable->GetX(0);
      G4double lastQ2 = theTable->GetX(tablePoints-1);
      G4double thePMax = 0;
      
      if (Qm2 > firstQ2)
	{
	  if (Qm2 < lastQ2)
	    {
	      //bisection to look for the index of Qm
	      size_t lowerBound = 0;
	      size_t upperBound = tablePoints-1;
	      while (lowerBound <= upperBound)
		{
		  size_t midBin = (lowerBound + upperBound)/2;
		  if( Qm2 < theTable->GetX(midBin))
		    { upperBound = midBin-1; }
		  else
		    { lowerBound = midBin+1; }
		}
	      //upperBound is the output (but also lowerBounf --> should be the same!)
	      G4double Q1 = theTable->GetX(upperBound);
 	      G4double Q2 = Qm2;
	      G4double DQ = (Q2-Q1)/((G4double)(nip-1));
	      G4double theA = theTable->GetA(upperBound);
	      G4double theB = theTable->GetB(upperBound);
	      G4double thePAC = theTable->GetPAC(upperBound);
	      G4DataVector* fun = new G4DataVector();
	      for (size_t k=0;k<nip;k++)
		{
		  G4double qi = Q1 + k*DQ;
		  G4double tau = (qi-Q1)/
		    (theTable->GetX(upperBound+1)-Q1);
		  G4double con1 = 2.0*theB*tau; 
		  G4double ci = 1.0+theA+theB;
		  G4double con2 = ci-theA*tau;
		  G4double etap = 0;
		  if (std::fabs(con1) > 1.0e-16*std::fabs(con2))
		    etap = con2*(1.0-std::sqrt(1.0-2.0*tau*con1/(con2*con2)))/con1;
		  else
		    etap = tau/con2;
		  G4double theFun = (theTable->GetPAC(upperBound+1)-thePAC)*
		    (1.0+(theA+theB*etap)*etap)*(1.0+(theA+theB*etap)*etap)/
		    ((1.0-theB*etap*etap)*ci*(theTable->GetX(upperBound+1)-Q1));
		  fun->push_back(theFun);
		}
	      //Now intergrate numerically the fun Cavalieri-Simpson's method
	      G4DataVector* sum = new G4DataVector;
	      G4double CONS = DQ*(1./12.);
	      G4double HCONS = 0.5*CONS;
	      sum->push_back(0.);
	      G4double secondPoint = (*sum)[0] + 
		(5.0*(*fun)[0]+8.0*(*fun)[1]-(*fun)[2])*CONS;
	      sum->push_back(secondPoint);
	      for (size_t hh=2;hh<nip-1;hh++)
		{
		  G4double previous = (*sum)[hh-1];
		  G4double next = previous+(13.0*((*fun)[hh-1]+(*fun)[hh])-
					    (*fun)[hh+1]-(*fun)[hh-2])*HCONS;
		  sum->push_back(next);
		}
	      G4double last = (*sum)[nip-2]+(5.0*(*fun)[nip-1]+8.0*(*fun)[nip-2]-
					     (*fun)[nip-3])*CONS;
	      sum->push_back(last);	 
	      thePMax = thePAC + (*sum)[sum->size()-1]; //last point
	      if (fun) delete fun;
	      if (sum) delete sum;
	    }
	  else
	    {
	      thePMax = 1.0;
	    }	 
	}
      else
	{
	  thePMax = theTable->GetPAC(0);
	}

      //Write number in the table
      theVec->PutValue(ie,energy,thePMax);
  }
  
  pMaxTable->insert(std::make_pair(mat,theVec));
  return;

}

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

void G4Penelope08RayleighModel::DumpFormFactorTable(const G4Material* mat)
{
  G4cout << "*****************************************************************" << G4endl;
  G4cout << "G4Penelope08RayleighModel: Form Factor Table for " << mat->GetName() << G4endl;
  //try to use the same format as Penelope-Fortran, namely Q (/m_e*c) and F
  G4cout <<  "Q/(m_e*c)                 F(Q)     " << G4endl;
  G4cout << "*****************************************************************" << G4endl;
  if (!logFormFactorTable->count(mat))
    BuildFormFactorTable(mat);
  
  G4PhysicsFreeVector* theVec = logFormFactorTable->find(mat)->second;
  for (size_t i=0;i<theVec->GetVectorLength();i++)
    {
      G4double logQ2 = theVec->GetLowEdgeEnergy(i);
      G4double Q = std::exp(0.5*logQ2);
      G4double logF2 = (*theVec)[i];
      G4double F = std::exp(0.5*logF2);
      G4cout << Q << "              " << F << G4endl;
    }
  //DONE
  return;
}