33 #define INCLXX_IN_GEANT4_MODE 1
46 namespace NuclearDensityFactory {
50 G4ThreadLocal std::map<G4int,NuclearDensity const *> *nuclearDensityCache = NULL;
51 G4ThreadLocal std::map<G4int,InverseInterpolationTable*> *rpCorrelationTableCache = NULL;
52 G4ThreadLocal std::map<G4int,InverseInterpolationTable*> *rCDFTableCache = NULL;
53 G4ThreadLocal std::map<G4int,InverseInterpolationTable*> *pCDFTableCache = NULL;
58 if(!nuclearDensityCache)
59 nuclearDensityCache =
new std::map<G4int,NuclearDensity const *>;
61 const G4int nuclideID = 1000*Z +
A;
62 const std::map<G4int,NuclearDensity const *>::const_iterator mapEntry = nuclearDensityCache->find(nuclideID);
63 if(mapEntry == nuclearDensityCache->end()) {
66 if(!rpCorrelationTableProton || !rpCorrelationTableNeutron)
69 (*nuclearDensityCache)[nuclideID] =
density;
72 return mapEntry->second;
79 if(!rpCorrelationTableCache)
80 rpCorrelationTableCache =
new std::map<G4int,InverseInterpolationTable*>;
82 const G4int nuclideID = ((t==
Proton) ? 1000 : -1000)*Z +
A;
83 const std::map<G4int,InverseInterpolationTable*>::const_iterator mapEntry = rpCorrelationTableCache->find(nuclideID);
84 if(mapEntry == rpCorrelationTableCache->end()) {
92 }
else if(A <= 19 && A > 6) {
97 }
else if(A <= 6 && A > 1) {
102 INCL_ERROR(
"No r-p correlation function for " << ((t==
Proton) ?
"protons" :
"neutrons") <<
" in A = "
103 << A <<
" Z = " << Z << std::endl);
113 normalisation(1./theFunction->integrate(xMin,xMax))
117 return Math::pow13(normalisation * theFunction->integrate(xMin,x));
122 } *theInverseCDFOneThird =
new InverseCDFOneThird(rpCorrelationFunction);
125 delete theInverseCDFOneThird;
126 delete rpCorrelationFunction;
127 INCL_DEBUG(
"Creating r-p correlation function for " << ((t==
Proton) ?
"protons" :
"neutrons") <<
" in A=" << A <<
", Z=" << Z <<
":"
128 << std::endl << theTable->
print() << std::endl);
130 (*rpCorrelationTableCache)[nuclideID] = theTable;
133 return mapEntry->second;
141 rCDFTableCache =
new std::map<G4int,InverseInterpolationTable*>;
143 const G4int nuclideID = ((t==
Proton) ? 1000 : -1000)*Z +
A;
144 const std::map<G4int,InverseInterpolationTable*>::const_iterator mapEntry = rCDFTableCache->find(nuclideID);
145 if(mapEntry == rCDFTableCache->end()) {
153 }
else if(A <= 19 && A > 6) {
158 }
else if(A <= 6 && A > 2) {
162 }
else if(A == 2 && Z == 1) {
165 INCL_ERROR(
"No nuclear density function for target A = "
166 << A <<
" Z = " << Z << std::endl);
171 delete rDensityFunction;
172 INCL_DEBUG(
"Creating inverse position CDF for A=" << A <<
", Z=" << Z <<
":" <<
173 std::endl << theTable->
print() << std::endl);
175 (*rCDFTableCache)[nuclideID] = theTable;
178 return mapEntry->second;
186 pCDFTableCache =
new std::map<G4int,InverseInterpolationTable*>;
188 const G4int nuclideID = ((t==
Proton) ? 1000 : -1000)*Z +
A;
189 const std::map<G4int,InverseInterpolationTable*>::const_iterator mapEntry = pCDFTableCache->find(nuclideID);
190 if(mapEntry == pCDFTableCache->end()) {
195 }
else if(A <= 19 && A > 2) {
198 }
else if(A == 2 && Z == 1) {
201 INCL_ERROR(
"No nuclear density function for target A = "
202 << A <<
" Z = " << Z << std::endl);
207 delete pDensityFunction;
208 INCL_DEBUG(
"Creating inverse momentum CDF for A=" << A <<
", Z=" << Z <<
":" <<
209 std::endl << theTable->
print() << std::endl);
211 (*pCDFTableCache)[nuclideID] = theTable;
214 return mapEntry->second;
221 if(!rpCorrelationTableCache)
222 rpCorrelationTableCache =
new std::map<G4int,InverseInterpolationTable*>;
224 const G4int nuclideID = ((t==
Proton) ? 1000 : -1000)*Z +
A;
225 const std::map<G4int,InverseInterpolationTable*>::const_iterator mapEntry = rpCorrelationTableCache->find(nuclideID);
226 if(mapEntry != rpCorrelationTableCache->end())
227 delete mapEntry->second;
229 (*rpCorrelationTableCache)[nuclideID] = table;
233 if(!nuclearDensityCache)
234 nuclearDensityCache =
new std::map<G4int,NuclearDensity const *>;
236 const G4int nuclideID = 1000*Z +
A;
237 const std::map<G4int,NuclearDensity const *>::const_iterator mapEntry = nuclearDensityCache->find(nuclideID);
238 if(mapEntry != nuclearDensityCache->end())
239 delete mapEntry->second;
241 (*nuclearDensityCache)[nuclideID] =
density;
246 if(nuclearDensityCache) {
247 for(std::map<G4int,NuclearDensity const *>::const_iterator i = nuclearDensityCache->begin(); i!=nuclearDensityCache->end(); ++i)
249 nuclearDensityCache->clear();
250 delete nuclearDensityCache;
251 nuclearDensityCache = NULL;
254 if(rpCorrelationTableCache) {
255 for(std::map<G4int,InverseInterpolationTable*>::const_iterator i = rpCorrelationTableCache->begin(); i!=rpCorrelationTableCache->end(); ++i)
257 rpCorrelationTableCache->clear();
258 delete rpCorrelationTableCache;
259 rpCorrelationTableCache = NULL;
263 for(std::map<G4int,InverseInterpolationTable*>::const_iterator i = rCDFTableCache->begin(); i!=rCDFTableCache->end(); ++i)
265 rCDFTableCache->clear();
266 delete rCDFTableCache;
267 rCDFTableCache = NULL;
271 for(std::map<G4int,InverseInterpolationTable*>::const_iterator i = pCDFTableCache->begin(); i!=pCDFTableCache->end(); ++i)
273 pCDFTableCache->clear();
274 delete pCDFTableCache;
275 pCDFTableCache = NULL;
Class for modified harmonic oscillator density.
InverseInterpolationTable * createPCDFTable(const ParticleType t, const G4int A, const G4int Z)
std::string print() const
NDF* class for the deuteron density according to the HardSphere potential.
InverseInterpolationTable * inverseCDFTable(const G4int nNodes=60) const
Return a pointer to the inverse of the CDF of this function.
Class for Gaussian density.
NuclearDensity const * createDensity(const G4int A, const G4int Z)
static const G4double A[nN]
const G4double oneOverSqrtThree
Class for interpolating the inverse of a 1-dimensional function.
G4double getSurfaceDiffuseness(const ParticleType t, const G4int A, const G4int Z)
G4double getRadiusParameter(const ParticleType t, const G4int A, const G4int Z)
G4double getMaximumNuclearRadius(const ParticleType t, const G4int A, const G4int Z)
void addDensityToCache(const G4int A, const G4int Z, NuclearDensity *const density)
NDF* class for the deuteron density according to the Paris potential.
InverseInterpolationTable * createRPCorrelationTable(const ParticleType t, const G4int A, const G4int Z)
InverseInterpolationTable * createRCDFTable(const ParticleType t, const G4int A, const G4int Z)
virtual G4double getXMinimum() const
Return the minimum allowed value of the independent variable.
G4double getMomentumRMS(const G4int A, const G4int Z)
Return the RMS of the momentum distribution (light clusters)
virtual G4double getXMaximum() const
Return the maximum allowed value of the independent variable.
G4double pow13(G4double x)
Class for Woods-Saxon density.
void addRPCorrelationToCache(const G4int A, const G4int Z, const ParticleType t, InverseInterpolationTable *const table)
G4ThreadLocal FermiMomentumFn getFermiMomentum