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()) {
69 if(!rpCorrelationTableProton || !rpCorrelationTableNeutron || !rpCorrelationTableLambda)
71 NuclearDensity const *density =
new NuclearDensity(A, Z, S, rpCorrelationTableProton, rpCorrelationTableNeutron, rpCorrelationTableLambda);
72 (*nuclearDensityCache)[nuclideID] = density;
75 return mapEntry->second;
82 if(!rpCorrelationTableCache)
83 rpCorrelationTableCache =
new std::map<G4int,InterpolationTable*>;
85 const G4int nuclideID = ((t==
Proton) ? 1000 : -1000)*Z +
A;
86 const std::map<G4int,InterpolationTable*>::const_iterator mapEntry = rpCorrelationTableCache->find(nuclideID);
87 if(mapEntry == rpCorrelationTableCache->end()) {
89 INCL_DEBUG(
"Creating r-p correlation function for " << ((t==
Proton) ?
"protons" :
"neutrons") <<
" in A=" << A <<
", Z=" << Z << std::endl);
97 INCL_DEBUG(
" ... Woods-Saxon; R0=" << radius <<
", a=" << diffuseness <<
", Rmax=" << maximumRadius << std::endl);
98 }
else if(A <= 19 && A > 6) {
103 INCL_DEBUG(
" ... MHO; param1=" << radius <<
", param2=" << diffuseness <<
", Rmax=" << maximumRadius << std::endl);
104 }
else if(A <= 6 && A > 1) {
108 INCL_DEBUG(
" ... Gaussian; sigma=" << radius <<
", Rmax=" << maximumRadius << std::endl);
110 INCL_ERROR(
"No r-p correlation function for " << ((t==
Proton) ?
"protons" :
"neutrons") <<
" in A = "
111 << A <<
" Z = " << Z <<
'\n');
116 delete rpCorrelationFunction;
117 INCL_DEBUG(
" ... here comes the table:\n" << theTable->
print() <<
'\n');
119 (*rpCorrelationTableCache)[nuclideID] = theTable;
122 return mapEntry->second;
130 rCDFTableCache =
new std::map<G4int,InterpolationTable*>;
132 const G4int nuclideID = ((t==
Proton) ? 1000 : -1000)*Z +
A;
133 const std::map<G4int,InterpolationTable*>::const_iterator mapEntry = rCDFTableCache->find(nuclideID);
134 if(mapEntry == rCDFTableCache->end()) {
142 }
else if(A <= 19 && A > 6) {
147 }
else if(A <= 6 && A > 2) {
151 }
else if(A == 2 && Z == 1) {
154 INCL_ERROR(
"No nuclear density function for target A = "
155 << A <<
" Z = " << Z <<
'\n');
160 delete rDensityFunction;
161 INCL_DEBUG(
"Creating inverse position CDF for A=" << A <<
", Z=" << Z <<
":" <<
162 '\n' << theTable->
print() <<
'\n');
164 (*rCDFTableCache)[nuclideID] = theTable;
167 return mapEntry->second;
175 pCDFTableCache =
new std::map<G4int,InterpolationTable*>;
177 const G4int nuclideID = ((t==
Proton) ? 1000 : -1000)*Z +
A;
178 const std::map<G4int,InterpolationTable*>::const_iterator mapEntry = pCDFTableCache->find(nuclideID);
179 if(mapEntry == pCDFTableCache->end()) {
184 }
else if(A <= 19 && A > 2) {
187 }
else if(A == 2 && Z == 1) {
190 INCL_ERROR(
"No nuclear density function for target A = "
191 << A <<
" Z = " << Z <<
'\n');
196 delete pDensityFunction;
197 INCL_DEBUG(
"Creating inverse momentum CDF for A=" << A <<
", Z=" << Z <<
":" <<
198 '\n' << theTable->
print() <<
'\n');
200 (*pCDFTableCache)[nuclideID] = theTable;
203 return mapEntry->second;
210 if(!rpCorrelationTableCache)
211 rpCorrelationTableCache =
new std::map<G4int,InterpolationTable*>;
213 const G4int nuclideID = ((t==
Proton) ? 1000 : -1000)*Z +
A;
214 const std::map<G4int,InterpolationTable*>::const_iterator mapEntry = rpCorrelationTableCache->find(nuclideID);
215 if(mapEntry != rpCorrelationTableCache->end())
216 delete mapEntry->second;
218 (*rpCorrelationTableCache)[nuclideID] = table;
222 if(!nuclearDensityCache)
223 nuclearDensityCache =
new std::map<G4int,NuclearDensity const *>;
225 const G4int nuclideID = 1000*Z +
A;
226 const std::map<G4int,NuclearDensity const *>::const_iterator mapEntry = nuclearDensityCache->find(nuclideID);
227 if(mapEntry != nuclearDensityCache->end())
228 delete mapEntry->second;
230 (*nuclearDensityCache)[nuclideID] = density;
235 if(nuclearDensityCache) {
236 for(std::map<G4int,NuclearDensity const *>::const_iterator i = nuclearDensityCache->begin(); i!=nuclearDensityCache->end(); ++i)
238 nuclearDensityCache->clear();
239 delete nuclearDensityCache;
240 nuclearDensityCache = NULL;
243 if(rpCorrelationTableCache) {
244 for(std::map<G4int,InterpolationTable*>::const_iterator i = rpCorrelationTableCache->begin(); i!=rpCorrelationTableCache->end(); ++i)
246 rpCorrelationTableCache->clear();
247 delete rpCorrelationTableCache;
248 rpCorrelationTableCache = NULL;
252 for(std::map<G4int,InterpolationTable*>::const_iterator i = rCDFTableCache->begin(); i!=rCDFTableCache->end(); ++i)
254 rCDFTableCache->clear();
255 delete rCDFTableCache;
256 rCDFTableCache = NULL;
260 for(std::map<G4int,InterpolationTable*>::const_iterator i = pCDFTableCache->begin(); i!=pCDFTableCache->end(); ++i)
262 pCDFTableCache->clear();
263 delete pCDFTableCache;
264 pCDFTableCache = NULL;