| Geant4 Cross Reference |
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 // INCL++ intra-nuclear cascade model
27 // Alain Boudard, CEA-Saclay, France
28 // Joseph Cugnon, University of Liege, Belgium
29 // Jean-Christophe David, CEA-Saclay, France
30 // Pekka Kaitaniemi, CEA-Saclay, France, and Helsinki Institute of Physics, Finland
31 // Sylvie Leray, CEA-Saclay, France
32 // Davide Mancusi, CEA-Saclay, France
33 //
34 #define INCLXX_IN_GEANT4_MODE 1
35
36 #include "globals.hh"
37
38 #include "G4INCLNuclearDensityFactory.hh"
39 #include "G4INCLNDFWoodsSaxon.hh"
40 #include "G4INCLNDFModifiedHarmonicOscillator.hh"
41 #include "G4INCLNDFGaussian.hh"
42 #include "G4INCLNDFParis.hh"
43 #include "G4INCLNDFHardSphere.hh"
44 #include "G4INCLInvFInterpolationTable.hh"
45
46 namespace G4INCL {
47
48 namespace NuclearDensityFactory {
49
50 namespace {
51
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;
56
57 }
58
59 NuclearDensity const *createDensity(const G4int A, const G4int Z, const G4int S) {
60 if(!nuclearDensityCache)
61 nuclearDensityCache = new std::map<G4int,NuclearDensity const *>;
62
63 const G4int nuclideID = 1000*Z + A; // MCNP-style nuclide IDs
64 const std::map<G4int,NuclearDensity const *>::const_iterator mapEntry = nuclearDensityCache->find(nuclideID);
65 if(mapEntry == nuclearDensityCache->end()) {
66 InterpolationTable *rpCorrelationTableProton = createRPCorrelationTable(Proton, A, Z);
67 InterpolationTable *rpCorrelationTableNeutron = createRPCorrelationTable(Neutron, A, Z);
68 InterpolationTable *rpCorrelationTableLambda = createRPCorrelationTable(Lambda, A, Z);
69 if(!rpCorrelationTableProton || !rpCorrelationTableNeutron || !rpCorrelationTableLambda)
70 return NULL;
71 NuclearDensity const *density = new NuclearDensity(A, Z, S, rpCorrelationTableProton, rpCorrelationTableNeutron, rpCorrelationTableLambda);
72 (*nuclearDensityCache)[nuclideID] = density;
73 return density;
74 } else {
75 return mapEntry->second;
76 }
77 }
78
79 InterpolationTable *createRPCorrelationTable(const ParticleType t, const G4int A, const G4int Z) {
80 // assert(t==Proton || t==Neutron || t==Lambda);
81
82 if(!rpCorrelationTableCache)
83 rpCorrelationTableCache = new std::map<G4int,InterpolationTable*>;
84
85 const G4int nuclideID = ((t==Proton) ? 1000 : -1000)*Z + A; // MCNP-style nuclide IDs
86 const std::map<G4int,InterpolationTable*>::const_iterator mapEntry = rpCorrelationTableCache->find(nuclideID);
87 if(mapEntry == rpCorrelationTableCache->end()) {
88
89 INCL_DEBUG("Creating r-p correlation function for " << ((t==Proton) ? "protons" : ((t==Neutron) ? "neutrons" : "lambdas")) << " in A=" << A << ", Z=" << Z << std::endl);
90
91 IFunction1D *rpCorrelationFunction;
92 if(A > 19) {
93 const G4double radius = ParticleTable::getRadiusParameter(t, A, Z);
94 const G4double diffuseness = ParticleTable::getSurfaceDiffuseness(t, A, Z);
95 const G4double maximumRadius = ParticleTable::getMaximumNuclearRadius(t, A, Z);
96 rpCorrelationFunction = new NuclearDensityFunctions::WoodsSaxonRP(radius, maximumRadius, diffuseness);
97 INCL_DEBUG(" ... Woods-Saxon; R0=" << radius << ", a=" << diffuseness << ", Rmax=" << maximumRadius << std::endl);
98 } else if(A <= 19 && A > 6) {
99 const G4double radius = ParticleTable::getRadiusParameter(t, A, Z);
100 const G4double diffuseness = ParticleTable::getSurfaceDiffuseness(t, A, Z);
101 const G4double maximumRadius = ParticleTable::getMaximumNuclearRadius(t, A, Z);
102 rpCorrelationFunction = new NuclearDensityFunctions::ModifiedHarmonicOscillatorRP(radius, maximumRadius, diffuseness);
103 INCL_DEBUG(" ... MHO; param1=" << radius << ", param2=" << diffuseness << ", Rmax=" << maximumRadius << std::endl);
104 } else if(A <= 6 && A > 1) { // Gaussian distribution for light nuclei
105 const G4double radius = ParticleTable::getRadiusParameter(t, A, Z);
106 const G4double maximumRadius = ParticleTable::getMaximumNuclearRadius(t, A, Z);
107 rpCorrelationFunction = new NuclearDensityFunctions::GaussianRP(maximumRadius, Math::oneOverSqrtThree * radius);
108 INCL_DEBUG(" ... Gaussian; sigma=" << radius << ", Rmax=" << maximumRadius << std::endl);
109 } else {
110 INCL_ERROR("No r-p correlation function for " << ((t==Proton) ? "protons" : "neutrons") << " in A = "
111 << A << " Z = " << Z << '\n');
112 return NULL;
113 }
114
115 InterpolationTable *theTable = rpCorrelationFunction->inverseCDFTable(Math::pow13);
116 delete rpCorrelationFunction;
117 INCL_DEBUG(" ... here comes the table:\n" << theTable->print() << '\n');
118
119 (*rpCorrelationTableCache)[nuclideID] = theTable;
120 return theTable;
121 } else {
122 return mapEntry->second;
123 }
124 }
125
126 InterpolationTable *createRCDFTable(const ParticleType t, const G4int A, const G4int Z) {
127 // assert(t==Proton || t==Neutron || t==Lambda);
128
129 if(!rCDFTableCache)
130 rCDFTableCache = new std::map<G4int,InterpolationTable*>;
131
132 const G4int nuclideID = ((t==Proton) ? 1000 : -1000)*Z + A; // MCNP-style nuclide IDs
133 const std::map<G4int,InterpolationTable*>::const_iterator mapEntry = rCDFTableCache->find(nuclideID);
134 if(mapEntry == rCDFTableCache->end()) {
135
136 IFunction1D *rDensityFunction;
137 if(A > 19) {
138 G4double radius = ParticleTable::getRadiusParameter(t, A, Z);
139 G4double diffuseness = ParticleTable::getSurfaceDiffuseness(t, A, Z);
140 G4double maximumRadius = ParticleTable::getMaximumNuclearRadius(t, A, Z);
141 rDensityFunction = new NuclearDensityFunctions::WoodsSaxon(radius, maximumRadius, diffuseness);
142 } else if(A <= 19 && A > 6) {
143 G4double radius = ParticleTable::getRadiusParameter(t, A, Z);
144 G4double diffuseness = ParticleTable::getSurfaceDiffuseness(t, A, Z);
145 G4double maximumRadius = ParticleTable::getMaximumNuclearRadius(t, A, Z);
146 rDensityFunction = new NuclearDensityFunctions::ModifiedHarmonicOscillator(radius, maximumRadius, diffuseness);
147 } else if(A <= 6 && A > 2) { // Gaussian distribution for light nuclei
148 G4double radius = ParticleTable::getRadiusParameter(t, A, Z);
149 G4double maximumRadius = ParticleTable::getMaximumNuclearRadius(t, A, Z);
150 rDensityFunction = new NuclearDensityFunctions::Gaussian(maximumRadius, Math::oneOverSqrtThree * radius);
151 } else if(A == 2 && Z == 1) { // density from the Paris potential for deuterons
152 rDensityFunction = new NuclearDensityFunctions::ParisR();
153 } else {
154 INCL_ERROR("No nuclear density function for target A = "
155 << A << " Z = " << Z << '\n');
156 return NULL;
157 }
158
159 InterpolationTable *theTable = rDensityFunction->inverseCDFTable();
160 delete rDensityFunction;
161 INCL_DEBUG("Creating inverse position CDF for A=" << A << ", Z=" << Z << ":" <<
162 '\n' << theTable->print() << '\n');
163
164 (*rCDFTableCache)[nuclideID] = theTable;
165 return theTable;
166 } else {
167 return mapEntry->second;
168 }
169 }
170
171 InterpolationTable *createPCDFTable(const ParticleType t, const G4int A, const G4int Z) {
172 // assert(t==Proton || t==Neutron || t==Lambda);
173
174 if(!pCDFTableCache)
175 pCDFTableCache = new std::map<G4int,InterpolationTable*>;
176
177 const G4int nuclideID = ((t==Proton) ? 1000 : -1000)*Z + A; // MCNP-style nuclide IDs
178 const std::map<G4int,InterpolationTable*>::const_iterator mapEntry = pCDFTableCache->find(nuclideID);
179 if(mapEntry == pCDFTableCache->end()) {
180 IFunction1D *pDensityFunction;
181 if(A > 19) {
182 const G4double theFermiMomentum = ParticleTable::getFermiMomentum(A, Z);
183 pDensityFunction = new NuclearDensityFunctions::HardSphere(theFermiMomentum);
184 } else if(A <= 19 && A > 2) { // Gaussian distribution for light nuclei
185 const G4double momentumRMS = Math::oneOverSqrtThree * ParticleTable::getMomentumRMS(A, Z);
186 pDensityFunction = new NuclearDensityFunctions::Gaussian(5.*momentumRMS, momentumRMS);
187 } else if(A == 2 && Z == 1) { // density from the Paris potential for deuterons
188 pDensityFunction = new NuclearDensityFunctions::ParisP();
189 } else {
190 INCL_ERROR("No nuclear density function for target A = "
191 << A << " Z = " << Z << '\n');
192 return NULL;
193 }
194
195 InterpolationTable *theTable = pDensityFunction->inverseCDFTable();
196 delete pDensityFunction;
197 INCL_DEBUG("Creating inverse momentum CDF for A=" << A << ", Z=" << Z << ":" <<
198 '\n' << theTable->print() << '\n');
199
200 (*pCDFTableCache)[nuclideID] = theTable;
201 return theTable;
202 } else {
203 return mapEntry->second;
204 }
205 }
206
207 void addRPCorrelationToCache(const G4int A, const G4int Z, const ParticleType t, InterpolationTable * const table) {
208 // assert(t==Proton || t==Neutron || t==Lambda);
209
210 if(!rpCorrelationTableCache)
211 rpCorrelationTableCache = new std::map<G4int,InterpolationTable*>;
212
213 const G4int nuclideID = ((t==Proton) ? 1000 : -1000)*Z + A; // MCNP-style nuclide IDs
214 const std::map<G4int,InterpolationTable*>::const_iterator mapEntry = rpCorrelationTableCache->find(nuclideID);
215 if(mapEntry != rpCorrelationTableCache->end())
216 delete mapEntry->second;
217
218 (*rpCorrelationTableCache)[nuclideID] = table;
219 }
220
221 void addDensityToCache(const G4int A, const G4int Z, NuclearDensity * const density) {
222 if(!nuclearDensityCache)
223 nuclearDensityCache = new std::map<G4int,NuclearDensity const *>;
224
225 const G4int nuclideID = 1000*Z + A; // MCNP-style nuclide IDs
226 const std::map<G4int,NuclearDensity const *>::const_iterator mapEntry = nuclearDensityCache->find(nuclideID);
227 if(mapEntry != nuclearDensityCache->end())
228 delete mapEntry->second;
229
230 (*nuclearDensityCache)[nuclideID] = density;
231 }
232
233 void clearCache() {
234
235 if(nuclearDensityCache) {
236 for(std::map<G4int,NuclearDensity const *>::const_iterator i = nuclearDensityCache->begin(); i!=nuclearDensityCache->end(); ++i)
237 delete i->second;
238 nuclearDensityCache->clear();
239 delete nuclearDensityCache;
240 nuclearDensityCache = NULL;
241 }
242
243 if(rpCorrelationTableCache) {
244 for(std::map<G4int,InterpolationTable*>::const_iterator i = rpCorrelationTableCache->begin(); i!=rpCorrelationTableCache->end(); ++i)
245 delete i->second;
246 rpCorrelationTableCache->clear();
247 delete rpCorrelationTableCache;
248 rpCorrelationTableCache = NULL;
249 }
250
251 if(rCDFTableCache) {
252 for(std::map<G4int,InterpolationTable*>::const_iterator i = rCDFTableCache->begin(); i!=rCDFTableCache->end(); ++i)
253 delete i->second;
254 rCDFTableCache->clear();
255 delete rCDFTableCache;
256 rCDFTableCache = NULL;
257 }
258
259 if(pCDFTableCache) {
260 for(std::map<G4int,InterpolationTable*>::const_iterator i = pCDFTableCache->begin(); i!=pCDFTableCache->end(); ++i)
261 delete i->second;
262 pCDFTableCache->clear();
263 delete pCDFTableCache;
264 pCDFTableCache = NULL;
265 }
266 }
267
268 } // namespace NuclearDensityFactory
269
270 } // namespace G4INCL
271