ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4INCLNuclearDensityFactory.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file 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, 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 
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" : "neutrons") << " in A=" << A << ", Z=" << Z << std::endl);
90 
91  IFunction1D *rpCorrelationFunction;
92  if(A > 19) {
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) {
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
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 
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) {
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) {
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
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 
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
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