lars/geant4
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity
Files
89a9605df1b6e52f0020a43d18e69f1293af21bf
geant4/source/processes/electromagnetic/lowenergy/src/G4Penelope08RayleighModel.cc
T
Gabriele Cosmo 89a9605df1 Import Geant4 9.4.0 source tree
2016-06-09 16:25:56 +02:00

1156 lines
37 KiB
C++
Raw Blame History

//
// ********************************************************************
// * 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 $
//
// 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;
}
Reference in New Issue View Git Blame Copy Permalink