G4INCL::NuclearDensityFactory Namespace Reference

Functions
InverseInterpolationTable *	createRPCorrelationTable (const ParticleType t, const G4int A, const G4int Z)
InverseInterpolationTable *	createRCDFTable (const ParticleType t, const G4int A, const G4int Z)
InverseInterpolationTable *	createPCDFTable (const ParticleType t, const G4int A, const G4int Z)
NuclearDensity const *	createDensity (const G4int A, const G4int Z)
void	addRPCorrelationToCache (const G4int A, const G4int Z, const ParticleType t, InverseInterpolationTable *const table)
void	addDensityToCache (const G4int A, const G4int Z, NuclearDensity *const density)
void	clearCache ()

Function Documentation

void G4INCL::NuclearDensityFactory::addDensityToCache	(	const G4int	A,
		const G4int	Z,
		NuclearDensity *const	density
	)

Definition at line 232 of file G4INCLNuclearDensityFactory.cc.

References density.

                                                                                         {
      if(!nuclearDensityCache)
        nuclearDensityCache = new std::map<G4int,NuclearDensity const *>;

      const G4int nuclideID = 1000*Z + A; // MCNP-style nuclide IDs
      const std::map<G4int,NuclearDensity const *>::const_iterator mapEntry = nuclearDensityCache->find(nuclideID);
      if(mapEntry != nuclearDensityCache->end())
        delete mapEntry->second;

      (*nuclearDensityCache)[nuclideID] = density;
    }

void G4INCL::NuclearDensityFactory::addRPCorrelationToCache	(	const G4int	A,
		const G4int	Z,
		const ParticleType	t,
		InverseInterpolationTable *const	table
	)

Definition at line 218 of file G4INCLNuclearDensityFactory.cc.

References G4INCL::Proton.

                                                                                                                              {
// assert(t==Proton || t==Neutron);

      if(!rpCorrelationTableCache)
        rpCorrelationTableCache = new std::map<G4int,InverseInterpolationTable*>;

      const G4int nuclideID = ((t==Proton) ? 1000 : -1000)*Z + A; // MCNP-style nuclide IDs
      const std::map<G4int,InverseInterpolationTable*>::const_iterator mapEntry = rpCorrelationTableCache->find(nuclideID);
      if(mapEntry != rpCorrelationTableCache->end())
        delete mapEntry->second;

      (*rpCorrelationTableCache)[nuclideID] = table;
    }

void G4INCL::NuclearDensityFactory::clearCache

(

)

Definition at line 244 of file G4INCLNuclearDensityFactory.cc.

Referenced by G4INCL::INCL::~INCL().

                      {

      if(nuclearDensityCache) {
        for(std::map<G4int,NuclearDensity const *>::const_iterator i = nuclearDensityCache->begin(); i!=nuclearDensityCache->end(); ++i)
          delete i->second;
        nuclearDensityCache->clear();
        delete nuclearDensityCache;
        nuclearDensityCache = NULL;
      }

      if(rpCorrelationTableCache) {
        for(std::map<G4int,InverseInterpolationTable*>::const_iterator i = rpCorrelationTableCache->begin(); i!=rpCorrelationTableCache->end(); ++i)
          delete i->second;
        rpCorrelationTableCache->clear();
        delete rpCorrelationTableCache;
        rpCorrelationTableCache = NULL;
      }

      if(rCDFTableCache) {
        for(std::map<G4int,InverseInterpolationTable*>::const_iterator i = rCDFTableCache->begin(); i!=rCDFTableCache->end(); ++i)
          delete i->second;
        rCDFTableCache->clear();
        delete rCDFTableCache;
        rCDFTableCache = NULL;
      }

      if(pCDFTableCache) {
        for(std::map<G4int,InverseInterpolationTable*>::const_iterator i = pCDFTableCache->begin(); i!=pCDFTableCache->end(); ++i)
          delete i->second;
        pCDFTableCache->clear();
        delete pCDFTableCache;
        pCDFTableCache = NULL;
      }
    }

NuclearDensity const * G4INCL::NuclearDensityFactory::createDensity	(	const G4int	A,
		const G4int	Z
	)

Definition at line 57 of file G4INCLNuclearDensityFactory.cc.

References createRPCorrelationTable(), density, G4INCL::Neutron, and G4INCL::Proton.

Referenced by G4INCL::Nucleus::Nucleus().

                                                                      {
      if(!nuclearDensityCache)
        nuclearDensityCache = new std::map<G4int,NuclearDensity const *>;

      const G4int nuclideID = 1000*Z + A; // MCNP-style nuclide IDs
      const std::map<G4int,NuclearDensity const *>::const_iterator mapEntry = nuclearDensityCache->find(nuclideID);
      if(mapEntry == nuclearDensityCache->end()) {
        InverseInterpolationTable *rpCorrelationTableProton = createRPCorrelationTable(Proton, A, Z);
        InverseInterpolationTable *rpCorrelationTableNeutron = createRPCorrelationTable(Neutron, A, Z);
        if(!rpCorrelationTableProton || !rpCorrelationTableNeutron)
          return NULL;
        NuclearDensity const *density = new NuclearDensity(A, Z, rpCorrelationTableProton, rpCorrelationTableNeutron);
        (*nuclearDensityCache)[nuclideID] = density;
        return density;
      } else {
        return mapEntry->second;
      }
    }

InverseInterpolationTable * G4INCL::NuclearDensityFactory::createPCDFTable	(	const ParticleType	t,
		const G4int	A,
		const G4int	Z
	)

Definition at line 182 of file G4INCLNuclearDensityFactory.cc.

References G4INCL::ParticleTable::getFermiMomentum, G4INCL::ParticleTable::getMomentumRMS(), INCL_DEBUG, INCL_ERROR, G4INCL::IFunction1D::inverseCDFTable(), G4INCL::Math::oneOverSqrtThree, G4INCL::InverseInterpolationTable::print(), and G4INCL::Proton.

Referenced by G4INCL::ParticleSampler::sampleParticles().

                                                                                                   {
// assert(t==Proton || t==Neutron);

      if(!pCDFTableCache)
        pCDFTableCache = new std::map<G4int,InverseInterpolationTable*>;

      const G4int nuclideID = ((t==Proton) ? 1000 : -1000)*Z + A; // MCNP-style nuclide IDs
      const std::map<G4int,InverseInterpolationTable*>::const_iterator mapEntry = pCDFTableCache->find(nuclideID);
      if(mapEntry == pCDFTableCache->end()) {
        IFunction1D *pDensityFunction;
        if(A > 19) {
          const G4double theFermiMomentum = ParticleTable::getFermiMomentum(A, Z);
          pDensityFunction = new NuclearDensityFunctions::HardSphere(theFermiMomentum);
        } else if(A <= 19 && A > 2) { // Gaussian distribution for light nuclei
          const G4double momentumRMS = Math::oneOverSqrtThree * ParticleTable::getMomentumRMS(A, Z);
          pDensityFunction = new NuclearDensityFunctions::Gaussian(5.*momentumRMS, momentumRMS);
        } else if(A == 2 && Z == 1) { // density from the Paris potential for deuterons
          pDensityFunction = new NuclearDensityFunctions::ParisP();
        } else {
          INCL_ERROR("No nuclear density function for target A = "
                << A << " Z = " << Z << std::endl);
          return NULL;
        }

        InverseInterpolationTable *theTable = pDensityFunction->inverseCDFTable();
        delete pDensityFunction;
        INCL_DEBUG("Creating inverse momentum CDF for A=" << A << ", Z=" << Z << ":" <<
              std::endl << theTable->print() << std::endl);

        (*pCDFTableCache)[nuclideID] = theTable;
        return theTable;
      } else {
        return mapEntry->second;
      }
    }

InverseInterpolationTable * G4INCL::NuclearDensityFactory::createRCDFTable	(	const ParticleType	t,
		const G4int	A,
		const G4int	Z
	)

Definition at line 137 of file G4INCLNuclearDensityFactory.cc.

References G4INCL::ParticleTable::getMaximumNuclearRadius(), G4INCL::ParticleTable::getRadiusParameter(), G4INCL::ParticleTable::getSurfaceDiffuseness(), INCL_DEBUG, INCL_ERROR, G4INCL::IFunction1D::inverseCDFTable(), G4INCL::Math::oneOverSqrtThree, G4INCL::InverseInterpolationTable::print(), and G4INCL::Proton.

Referenced by G4INCL::ParticleSampler::sampleParticles().

                                                                                                   {
// assert(t==Proton || t==Neutron);

      if(!rCDFTableCache)
        rCDFTableCache = new std::map<G4int,InverseInterpolationTable*>;

      const G4int nuclideID = ((t==Proton) ? 1000 : -1000)*Z + A; // MCNP-style nuclide IDs
      const std::map<G4int,InverseInterpolationTable*>::const_iterator mapEntry = rCDFTableCache->find(nuclideID);
      if(mapEntry == rCDFTableCache->end()) {

        IFunction1D *rDensityFunction;
        if(A > 19) {
          G4double radius = ParticleTable::getRadiusParameter(t, A, Z);
          G4double diffuseness = ParticleTable::getSurfaceDiffuseness(t, A, Z);
          G4double maximumRadius = ParticleTable::getMaximumNuclearRadius(t, A, Z);
          rDensityFunction = new NuclearDensityFunctions::WoodsSaxon(radius, maximumRadius, diffuseness);
        } else if(A <= 19 && A > 6) {
          G4double radius = ParticleTable::getRadiusParameter(t, A, Z);
          G4double diffuseness = ParticleTable::getSurfaceDiffuseness(t, A, Z);
          G4double maximumRadius = ParticleTable::getMaximumNuclearRadius(t, A, Z);
          rDensityFunction = new NuclearDensityFunctions::ModifiedHarmonicOscillator(radius, maximumRadius, diffuseness);
        } else if(A <= 6 && A > 2) { // Gaussian distribution for light nuclei
          G4double radius = ParticleTable::getRadiusParameter(t, A, Z);
          G4double maximumRadius = ParticleTable::getMaximumNuclearRadius(t, A, Z);
          rDensityFunction = new NuclearDensityFunctions::Gaussian(maximumRadius, Math::oneOverSqrtThree * radius);
        } else if(A == 2 && Z == 1) { // density from the Paris potential for deuterons
          rDensityFunction = new NuclearDensityFunctions::ParisR();
        } else {
          INCL_ERROR("No nuclear density function for target A = "
                << A << " Z = " << Z << std::endl);
          return NULL;
        }

        InverseInterpolationTable *theTable = rDensityFunction->inverseCDFTable();
        delete rDensityFunction;
        INCL_DEBUG("Creating inverse position CDF for A=" << A << ", Z=" << Z << ":" <<
              std::endl << theTable->print() << std::endl);

        (*rCDFTableCache)[nuclideID] = theTable;
        return theTable;
      } else {
        return mapEntry->second;
      }
    }

InverseInterpolationTable * G4INCL::NuclearDensityFactory::createRPCorrelationTable	(	const ParticleType	t,
		const G4int	A,
		const G4int	Z
	)

Definition at line 76 of file G4INCLNuclearDensityFactory.cc.

References G4INCL::ParticleTable::getMaximumNuclearRadius(), G4INCL::ParticleTable::getRadiusParameter(), G4INCL::ParticleTable::getSurfaceDiffuseness(), G4INCL::IFunction1D::getXMaximum(), G4INCL::IFunction1D::getXMinimum(), INCL_DEBUG, INCL_ERROR, G4INCL::Math::oneOverSqrtThree, G4INCL::Math::pow13(), G4INCL::InverseInterpolationTable::print(), G4INCL::Proton, and test::x.

Referenced by createDensity().

                                                                                                            {
// assert(t==Proton || t==Neutron);

      if(!rpCorrelationTableCache)
        rpCorrelationTableCache = new std::map<G4int,InverseInterpolationTable*>;

      const G4int nuclideID = ((t==Proton) ? 1000 : -1000)*Z + A; // MCNP-style nuclide IDs
      const std::map<G4int,InverseInterpolationTable*>::const_iterator mapEntry = rpCorrelationTableCache->find(nuclideID);
      if(mapEntry == rpCorrelationTableCache->end()) {

        IFunction1D *rpCorrelationFunction;
        if(A > 19) {
          G4double radius = ParticleTable::getRadiusParameter(t, A, Z);
          G4double diffuseness = ParticleTable::getSurfaceDiffuseness(t, A, Z);
          const G4double maximumRadius = ParticleTable::getMaximumNuclearRadius(t, A, Z);
          rpCorrelationFunction = new NuclearDensityFunctions::WoodsSaxonRP(radius, maximumRadius, diffuseness);
        } else if(A <= 19 && A > 6) {
          const G4double radius = ParticleTable::getRadiusParameter(t, A, Z);
          const G4double diffuseness = ParticleTable::getSurfaceDiffuseness(t, A, Z);
          const G4double maximumRadius = ParticleTable::getMaximumNuclearRadius(t, A, Z);
          rpCorrelationFunction = new NuclearDensityFunctions::ModifiedHarmonicOscillatorRP(radius, maximumRadius, diffuseness);
        } else if(A <= 6 && A > 1) { // Gaussian distribution for light nuclei
          const G4double radius = ParticleTable::getRadiusParameter(t, A, Z);
          const G4double maximumRadius = ParticleTable::getMaximumNuclearRadius(t, A, Z);
          rpCorrelationFunction = new NuclearDensityFunctions::GaussianRP(maximumRadius, Math::oneOverSqrtThree * radius);
        } else {
          INCL_ERROR("No r-p correlation function for " << ((t==Proton) ? "protons" : "neutrons") << " in A = "
                << A << " Z = " << Z << std::endl);
          return NULL;

        }

        class InverseCDFOneThird : public IFunction1D {
          public:
            InverseCDFOneThird(IFunction1D const * const f) :
              IFunction1D(f->getXMinimum(), f->getXMaximum()),
              theFunction(f),
              normalisation(1./theFunction->integrate(xMin,xMax))
          {}

            G4double operator()(const G4double x) const {
              return Math::pow13(normalisation * theFunction->integrate(xMin,x));
            }
          private:
            IFunction1D const * const theFunction;
            const G4double normalisation;
        } *theInverseCDFOneThird = new InverseCDFOneThird(rpCorrelationFunction);

        InverseInterpolationTable *theTable = new InverseInterpolationTable(*theInverseCDFOneThird);
        delete theInverseCDFOneThird;
        delete rpCorrelationFunction;
        INCL_DEBUG("Creating r-p correlation function for " << ((t==Proton) ? "protons" : "neutrons") << " in A=" << A << ", Z=" << Z << ":"
              << std::endl << theTable->print() << std::endl);

        (*rpCorrelationTableCache)[nuclideID] = theTable;
        return theTable;
      } else {
        return mapEntry->second;
      }
    }