Geant4  v4-10.4-release
 모두 클래스 네임스페이스들 파일들 함수 변수 타입정의 열거형 타입 열거형 멤버 Friends 매크로 그룹들 페이지들
G4INCLNuclearDensityFactory.cc
이 파일의 문서화 페이지로 가기
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 
39 #include "G4INCLNDFWoodsSaxon.hh"
41 #include "G4INCLNDFGaussian.hh"
42 #include "G4INCLNDFParis.hh"
43 #include "G4INCLNDFHardSphere.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) {
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  if(!rpCorrelationTableProton || !rpCorrelationTableNeutron)
69  return NULL;
70  NuclearDensity const *density = new NuclearDensity(A, Z, rpCorrelationTableProton, rpCorrelationTableNeutron);
71  (*nuclearDensityCache)[nuclideID] = density;
72  return density;
73  } else {
74  return mapEntry->second;
75  }
76  }
77 
79 // assert(t==Proton || t==Neutron);
80 
81  if(!rpCorrelationTableCache)
82  rpCorrelationTableCache = new std::map<G4int,InterpolationTable*>;
83 
84  const G4int nuclideID = ((t==Proton) ? 1000 : -1000)*Z + A; // MCNP-style nuclide IDs
85  const std::map<G4int,InterpolationTable*>::const_iterator mapEntry = rpCorrelationTableCache->find(nuclideID);
86  if(mapEntry == rpCorrelationTableCache->end()) {
87 
88  INCL_DEBUG("Creating r-p correlation function for " << ((t==Proton) ? "protons" : "neutrons") << " in A=" << A << ", Z=" << Z << std::endl);
89 
90  IFunction1D *rpCorrelationFunction;
91  if(A > 19) {
93  const G4double diffuseness = ParticleTable::getSurfaceDiffuseness(t, A, Z);
94  const G4double maximumRadius = ParticleTable::getMaximumNuclearRadius(t, A, Z);
95  rpCorrelationFunction = new NuclearDensityFunctions::WoodsSaxonRP(radius, maximumRadius, diffuseness);
96  INCL_DEBUG(" ... Woods-Saxon; R0=" << radius << ", a=" << diffuseness << ", Rmax=" << maximumRadius << std::endl);
97  } else if(A <= 19 && A > 6) {
99  const G4double diffuseness = ParticleTable::getSurfaceDiffuseness(t, A, Z);
100  const G4double maximumRadius = ParticleTable::getMaximumNuclearRadius(t, A, Z);
101  rpCorrelationFunction = new NuclearDensityFunctions::ModifiedHarmonicOscillatorRP(radius, maximumRadius, diffuseness);
102  INCL_DEBUG(" ... MHO; param1=" << radius << ", param2=" << diffuseness << ", Rmax=" << maximumRadius << std::endl);
103  } else if(A <= 6 && A > 1) { // Gaussian distribution for light nuclei
105  const G4double maximumRadius = ParticleTable::getMaximumNuclearRadius(t, A, Z);
106  rpCorrelationFunction = new NuclearDensityFunctions::GaussianRP(maximumRadius, Math::oneOverSqrtThree * radius);
107  INCL_DEBUG(" ... Gaussian; sigma=" << radius << ", Rmax=" << maximumRadius << std::endl);
108  } else {
109  INCL_ERROR("No r-p correlation function for " << ((t==Proton) ? "protons" : "neutrons") << " in A = "
110  << A << " Z = " << Z << '\n');
111  return NULL;
112  }
113 
114  InterpolationTable *theTable = rpCorrelationFunction->inverseCDFTable(Math::pow13);
115  delete rpCorrelationFunction;
116  INCL_DEBUG(" ... here comes the table:\n" << theTable->print() << '\n');
117 
118  (*rpCorrelationTableCache)[nuclideID] = theTable;
119  return theTable;
120  } else {
121  return mapEntry->second;
122  }
123  }
124 
126 // assert(t==Proton || t==Neutron);
127 
128  if(!rCDFTableCache)
129  rCDFTableCache = new std::map<G4int,InterpolationTable*>;
130 
131  const G4int nuclideID = ((t==Proton) ? 1000 : -1000)*Z + A; // MCNP-style nuclide IDs
132  const std::map<G4int,InterpolationTable*>::const_iterator mapEntry = rCDFTableCache->find(nuclideID);
133  if(mapEntry == rCDFTableCache->end()) {
134 
135  IFunction1D *rDensityFunction;
136  if(A > 19) {
138  G4double diffuseness = ParticleTable::getSurfaceDiffuseness(t, A, Z);
139  G4double maximumRadius = ParticleTable::getMaximumNuclearRadius(t, A, Z);
140  rDensityFunction = new NuclearDensityFunctions::WoodsSaxon(radius, maximumRadius, diffuseness);
141  } else if(A <= 19 && A > 6) {
143  G4double diffuseness = ParticleTable::getSurfaceDiffuseness(t, A, Z);
144  G4double maximumRadius = ParticleTable::getMaximumNuclearRadius(t, A, Z);
145  rDensityFunction = new NuclearDensityFunctions::ModifiedHarmonicOscillator(radius, maximumRadius, diffuseness);
146  } else if(A <= 6 && A > 2) { // Gaussian distribution for light nuclei
148  G4double maximumRadius = ParticleTable::getMaximumNuclearRadius(t, A, Z);
149  rDensityFunction = new NuclearDensityFunctions::Gaussian(maximumRadius, Math::oneOverSqrtThree * radius);
150  } else if(A == 2 && Z == 1) { // density from the Paris potential for deuterons
151  rDensityFunction = new NuclearDensityFunctions::ParisR();
152  } else {
153  INCL_ERROR("No nuclear density function for target A = "
154  << A << " Z = " << Z << '\n');
155  return NULL;
156  }
157 
158  InterpolationTable *theTable = rDensityFunction->inverseCDFTable();
159  delete rDensityFunction;
160  INCL_DEBUG("Creating inverse position CDF for A=" << A << ", Z=" << Z << ":" <<
161  '\n' << theTable->print() << '\n');
162 
163  (*rCDFTableCache)[nuclideID] = theTable;
164  return theTable;
165  } else {
166  return mapEntry->second;
167  }
168  }
169 
171 // assert(t==Proton || t==Neutron);
172 
173  if(!pCDFTableCache)
174  pCDFTableCache = new std::map<G4int,InterpolationTable*>;
175 
176  const G4int nuclideID = ((t==Proton) ? 1000 : -1000)*Z + A; // MCNP-style nuclide IDs
177  const std::map<G4int,InterpolationTable*>::const_iterator mapEntry = pCDFTableCache->find(nuclideID);
178  if(mapEntry == pCDFTableCache->end()) {
179  IFunction1D *pDensityFunction;
180  if(A > 19) {
181  const G4double theFermiMomentum = ParticleTable::getFermiMomentum(A, Z);
182  pDensityFunction = new NuclearDensityFunctions::HardSphere(theFermiMomentum);
183  } else if(A <= 19 && A > 2) { // Gaussian distribution for light nuclei
185  pDensityFunction = new NuclearDensityFunctions::Gaussian(5.*momentumRMS, momentumRMS);
186  } else if(A == 2 && Z == 1) { // density from the Paris potential for deuterons
187  pDensityFunction = new NuclearDensityFunctions::ParisP();
188  } else {
189  INCL_ERROR("No nuclear density function for target A = "
190  << A << " Z = " << Z << '\n');
191  return NULL;
192  }
193 
194  InterpolationTable *theTable = pDensityFunction->inverseCDFTable();
195  delete pDensityFunction;
196  INCL_DEBUG("Creating inverse momentum CDF for A=" << A << ", Z=" << Z << ":" <<
197  '\n' << theTable->print() << '\n');
198 
199  (*pCDFTableCache)[nuclideID] = theTable;
200  return theTable;
201  } else {
202  return mapEntry->second;
203  }
204  }
205 
206  void addRPCorrelationToCache(const G4int A, const G4int Z, const ParticleType t, InterpolationTable * const table) {
207 // assert(t==Proton || t==Neutron);
208 
209  if(!rpCorrelationTableCache)
210  rpCorrelationTableCache = new std::map<G4int,InterpolationTable*>;
211 
212  const G4int nuclideID = ((t==Proton) ? 1000 : -1000)*Z + A; // MCNP-style nuclide IDs
213  const std::map<G4int,InterpolationTable*>::const_iterator mapEntry = rpCorrelationTableCache->find(nuclideID);
214  if(mapEntry != rpCorrelationTableCache->end())
215  delete mapEntry->second;
216 
217  (*rpCorrelationTableCache)[nuclideID] = table;
218  }
219 
220  void addDensityToCache(const G4int A, const G4int Z, NuclearDensity * const density) {
221  if(!nuclearDensityCache)
222  nuclearDensityCache = new std::map<G4int,NuclearDensity const *>;
223 
224  const G4int nuclideID = 1000*Z + A; // MCNP-style nuclide IDs
225  const std::map<G4int,NuclearDensity const *>::const_iterator mapEntry = nuclearDensityCache->find(nuclideID);
226  if(mapEntry != nuclearDensityCache->end())
227  delete mapEntry->second;
228 
229  (*nuclearDensityCache)[nuclideID] = density;
230  }
231 
232  void clearCache() {
233 
234  if(nuclearDensityCache) {
235  for(std::map<G4int,NuclearDensity const *>::const_iterator i = nuclearDensityCache->begin(); i!=nuclearDensityCache->end(); ++i)
236  delete i->second;
237  nuclearDensityCache->clear();
238  delete nuclearDensityCache;
239  nuclearDensityCache = NULL;
240  }
241 
242  if(rpCorrelationTableCache) {
243  for(std::map<G4int,InterpolationTable*>::const_iterator i = rpCorrelationTableCache->begin(); i!=rpCorrelationTableCache->end(); ++i)
244  delete i->second;
245  rpCorrelationTableCache->clear();
246  delete rpCorrelationTableCache;
247  rpCorrelationTableCache = NULL;
248  }
249 
250  if(rCDFTableCache) {
251  for(std::map<G4int,InterpolationTable*>::const_iterator i = rCDFTableCache->begin(); i!=rCDFTableCache->end(); ++i)
252  delete i->second;
253  rCDFTableCache->clear();
254  delete rCDFTableCache;
255  rCDFTableCache = NULL;
256  }
257 
258  if(pCDFTableCache) {
259  for(std::map<G4int,InterpolationTable*>::const_iterator i = pCDFTableCache->begin(); i!=pCDFTableCache->end(); ++i)
260  delete i->second;
261  pCDFTableCache->clear();
262  delete pCDFTableCache;
263  pCDFTableCache = NULL;
264  }
265  }
266 
267  } // namespace NuclearDensityFactory
268 
269 } // namespace G4INCL
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)
#define G4ThreadLocal
Definition: tls.hh:69
InterpolationTable * createRCDFTable(const ParticleType t, const G4int A, const G4int Z)
#define INCL_DEBUG(x)
Float_t Z
G4double getMaximumNuclearRadius(const ParticleType t, const G4int A, const G4int Z)
double G4double
Definition: G4Types.hh:76
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.
#define INCL_ERROR(x)
G4double pow13(G4double x)
double A(double temperature)
Double_t radius
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)
int G4int
Definition: G4Types.hh:78
NDF* class for the deuteron density according to the Paris potential.
const G4double oneOverSqrtThree
Class for Woods-Saxon density.