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 * createPCDFTable(const ParticleType t, const G4int A, const G4int Z)
Class for modified harmonic oscillator density.
Simple interpolation table for the inverse of a IFunction1D functor.
void addRPCorrelationToCache(const G4int A, const G4int Z, const ParticleType t, InterpolationTable *const table)
InterpolationTable * createRCDFTable(const ParticleType t, const G4int A, const G4int Z)
NDF* class for the deuteron density according to the HardSphere potential.
std::string print() const
double A(double temperature)
InterpolationTable * createRPCorrelationTable(const ParticleType t, const G4int A, const G4int Z)
Class for Gaussian density.
NuclearDensity const * createDensity(const G4int A, const G4int Z)
const G4double oneOverSqrtThree
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.
G4double getMomentumRMS(const G4int A, const G4int Z)
Return the RMS of the momentum distribution (light clusters)
G4double pow13(G4double x)
InterpolationTable * inverseCDFTable(ManipulatorFunc fWrap=0, const G4int nNodes=60) const
Return a pointer to the inverse of the CDF of this function.
Class for Woods-Saxon density.
Class for interpolating the of a 1-dimensional function.
G4ThreadLocal FermiMomentumFn getFermiMomentum