34 #define INCLXX_IN_GEANT4_MODE 1
48 namespace NuclearDensityFactory {
52 G4ThreadLocal std::map<G4int,NuclearDensity const *> *nuclearDensityCache = NULL;
53 G4ThreadLocal std::map<G4int,InterpolationTable*> *rpCorrelationTableCache = NULL;
54 G4ThreadLocal std::map<G4int,InterpolationTable*> *rCDFTableCache = NULL;
55 G4ThreadLocal std::map<G4int,InterpolationTable*> *pCDFTableCache = NULL;
60 if(!nuclearDensityCache)
61 nuclearDensityCache =
new std::map<G4int,NuclearDensity const *>;
63 const G4int nuclideID = 1000*Z +
A;
64 const std::map<G4int,NuclearDensity const *>::const_iterator mapEntry = nuclearDensityCache->find(nuclideID);
65 if(mapEntry == nuclearDensityCache->end()) {
68 if(!rpCorrelationTableProton || !rpCorrelationTableNeutron)
71 (*nuclearDensityCache)[nuclideID] = density;
74 return mapEntry->second;
81 if(!rpCorrelationTableCache)
82 rpCorrelationTableCache =
new std::map<G4int,InterpolationTable*>;
84 const G4int nuclideID = ((t==
Proton) ? 1000 : -1000)*Z +
A;
85 const std::map<G4int,InterpolationTable*>::const_iterator mapEntry = rpCorrelationTableCache->find(nuclideID);
86 if(mapEntry == rpCorrelationTableCache->end()) {
88 INCL_DEBUG(
"Creating r-p correlation function for " << ((t==
Proton) ?
"protons" :
"neutrons") <<
" in A=" << A <<
", Z=" << Z << std::endl);
96 INCL_DEBUG(
" ... Woods-Saxon; R0=" << radius <<
", a=" << diffuseness <<
", Rmax=" << maximumRadius << std::endl);
97 }
else if(A <= 19 && A > 6) {
102 INCL_DEBUG(
" ... MHO; param1=" << radius <<
", param2=" << diffuseness <<
", Rmax=" << maximumRadius << std::endl);
103 }
else if(A <= 6 && A > 1) {
107 INCL_DEBUG(
" ... Gaussian; sigma=" << radius <<
", Rmax=" << maximumRadius << std::endl);
109 INCL_ERROR(
"No r-p correlation function for " << ((t==
Proton) ?
"protons" :
"neutrons") <<
" in A = "
110 << A <<
" Z = " << Z <<
'\n');
115 delete rpCorrelationFunction;
116 INCL_DEBUG(
" ... here comes the table:\n" << theTable->
print() <<
'\n');
118 (*rpCorrelationTableCache)[nuclideID] = theTable;
121 return mapEntry->second;
129 rCDFTableCache =
new std::map<G4int,InterpolationTable*>;
131 const G4int nuclideID = ((t==
Proton) ? 1000 : -1000)*Z +
A;
132 const std::map<G4int,InterpolationTable*>::const_iterator mapEntry = rCDFTableCache->find(nuclideID);
133 if(mapEntry == rCDFTableCache->end()) {
141 }
else if(A <= 19 && A > 6) {
146 }
else if(A <= 6 && A > 2) {
150 }
else if(A == 2 && Z == 1) {
153 INCL_ERROR(
"No nuclear density function for target A = "
154 << A <<
" Z = " << Z <<
'\n');
159 delete rDensityFunction;
160 INCL_DEBUG(
"Creating inverse position CDF for A=" << A <<
", Z=" << Z <<
":" <<
161 '\n' << theTable->
print() <<
'\n');
163 (*rCDFTableCache)[nuclideID] = theTable;
166 return mapEntry->second;
174 pCDFTableCache =
new std::map<G4int,InterpolationTable*>;
176 const G4int nuclideID = ((t==
Proton) ? 1000 : -1000)*Z +
A;
177 const std::map<G4int,InterpolationTable*>::const_iterator mapEntry = pCDFTableCache->find(nuclideID);
178 if(mapEntry == pCDFTableCache->end()) {
183 }
else if(A <= 19 && A > 2) {
186 }
else if(A == 2 && Z == 1) {
189 INCL_ERROR(
"No nuclear density function for target A = "
190 << A <<
" Z = " << Z <<
'\n');
195 delete pDensityFunction;
196 INCL_DEBUG(
"Creating inverse momentum CDF for A=" << A <<
", Z=" << Z <<
":" <<
197 '\n' << theTable->
print() <<
'\n');
199 (*pCDFTableCache)[nuclideID] = theTable;
202 return mapEntry->second;
209 if(!rpCorrelationTableCache)
210 rpCorrelationTableCache =
new std::map<G4int,InterpolationTable*>;
212 const G4int nuclideID = ((t==
Proton) ? 1000 : -1000)*Z +
A;
213 const std::map<G4int,InterpolationTable*>::const_iterator mapEntry = rpCorrelationTableCache->find(nuclideID);
214 if(mapEntry != rpCorrelationTableCache->end())
215 delete mapEntry->second;
217 (*rpCorrelationTableCache)[nuclideID] = table;
221 if(!nuclearDensityCache)
222 nuclearDensityCache =
new std::map<G4int,NuclearDensity const *>;
224 const G4int nuclideID = 1000*Z +
A;
225 const std::map<G4int,NuclearDensity const *>::const_iterator mapEntry = nuclearDensityCache->find(nuclideID);
226 if(mapEntry != nuclearDensityCache->end())
227 delete mapEntry->second;
229 (*nuclearDensityCache)[nuclideID] = density;
234 if(nuclearDensityCache) {
235 for(std::map<G4int,NuclearDensity const *>::const_iterator i = nuclearDensityCache->begin(); i!=nuclearDensityCache->end(); ++i)
237 nuclearDensityCache->clear();
238 delete nuclearDensityCache;
239 nuclearDensityCache = NULL;
242 if(rpCorrelationTableCache) {
243 for(std::map<G4int,InterpolationTable*>::const_iterator i = rpCorrelationTableCache->begin(); i!=rpCorrelationTableCache->end(); ++i)
245 rpCorrelationTableCache->clear();
246 delete rpCorrelationTableCache;
247 rpCorrelationTableCache = NULL;
251 for(std::map<G4int,InterpolationTable*>::const_iterator i = rCDFTableCache->begin(); i!=rCDFTableCache->end(); ++i)
253 rCDFTableCache->clear();
254 delete rCDFTableCache;
255 rCDFTableCache = NULL;
259 for(std::map<G4int,InterpolationTable*>::const_iterator i = pCDFTableCache->begin(); i!=pCDFTableCache->end(); ++i)
261 pCDFTableCache->clear();
262 delete pCDFTableCache;
263 pCDFTableCache = NULL;
InterpolationTable * inverseCDFTable(ManipulatorFunc fWrap=0, const G4int nNodes=60) const
Return a pointer to the inverse of the CDF of this function.
G4double getRadiusParameter(const ParticleType t, const G4int A, const G4int Z)
void addRPCorrelationToCache(const G4int A, const G4int Z, const ParticleType t, InterpolationTable *const table)
G4double getMomentumRMS(const G4int A, const G4int Z)
Return the RMS of the momentum distribution (light clusters)
G4ThreadLocal FermiMomentumFn getFermiMomentum
InterpolationTable * createRPCorrelationTable(const ParticleType t, const G4int A, const G4int Z)
NuclearDensity const * createDensity(const G4int A, const G4int Z)
InterpolationTable * createRCDFTable(const ParticleType t, const G4int A, const G4int Z)
G4double getMaximumNuclearRadius(const ParticleType t, const G4int A, const G4int Z)
NDF* class for the deuteron density according to the HardSphere potential.
Class for modified harmonic oscillator density.
G4double getSurfaceDiffuseness(const ParticleType t, const G4int A, const G4int Z)
Class for Gaussian density.
G4double pow13(G4double x)
double A(double temperature)
void addDensityToCache(const G4int A, const G4int Z, NuclearDensity *const density)
Class for interpolating the of a 1-dimensional function.
Simple interpolation table for the inverse of a IFunction1D functor.
InterpolationTable * createPCDFTable(const ParticleType t, const G4int A, const G4int Z)
NDF* class for the deuteron density according to the Paris potential.
std::string print() const
const G4double oneOverSqrtThree
Class for Woods-Saxon density.