00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033 #define INCLXX_IN_GEANT4_MODE 1
00034
00035 #include "globals.hh"
00036
00037 #include "G4INCLNuclearDensityFactory.hh"
00038 #include "G4INCLNDFWoodsSaxon.hh"
00039 #include "G4INCLNDFModifiedHarmonicOscillator.hh"
00040 #include "G4INCLNDFGaussian.hh"
00041 #include "G4INCLNDFParis.hh"
00042 #include "G4INCLNDFHardSphere.hh"
00043
00044 namespace G4INCL {
00045
00046 std::map<G4int,NuclearDensity*> NuclearDensityFactory::nuclearDensityCache;
00047
00048 std::map<G4int,InverseInterpolationTable*> NuclearDensityFactory::rpCorrelationTableCache;
00049
00050 std::map<G4int,InverseInterpolationTable*> NuclearDensityFactory::rCDFTableCache;
00051
00052 std::map<G4int,InverseInterpolationTable*> NuclearDensityFactory::pCDFTableCache;
00053
00054 NuclearDensity* NuclearDensityFactory::createDensity(const G4int A, const G4int Z) {
00055 const G4int nuclideID = 1000*Z + A;
00056 const std::map<G4int,NuclearDensity*>::const_iterator mapEntry = nuclearDensityCache.find(nuclideID);
00057 if(mapEntry == nuclearDensityCache.end()) {
00058 InverseInterpolationTable *rpCorrelationTable = NuclearDensityFactory::createRPCorrelationTable(A, Z);
00059 if(!rpCorrelationTable)
00060 return NULL;
00061 NuclearDensity *density = new NuclearDensity(A, Z, rpCorrelationTable);
00062 nuclearDensityCache[nuclideID] = density;
00063 return density;
00064 } else {
00065 return mapEntry->second;
00066 }
00067 }
00068
00069 InverseInterpolationTable *NuclearDensityFactory::createRPCorrelationTable(const G4int A, const G4int Z) {
00070 const G4int nuclideID = 1000*Z + A;
00071 const std::map<G4int,InverseInterpolationTable*>::const_iterator mapEntry = rpCorrelationTableCache.find(nuclideID);
00072 if(mapEntry == rpCorrelationTableCache.end()) {
00073
00074 IFunction1D *rpCorrelationFunction;
00075 if(A > 19) {
00076 const G4double radius = ParticleTable::getRadiusParameter(A, Z);
00077 const G4double diffuseness = ParticleTable::getSurfaceDiffuseness(A, Z);
00078 const G4double maximumRadius = ParticleTable::getMaximumNuclearRadius(A, Z);
00079 rpCorrelationFunction = new NuclearDensityFunctions::WoodsSaxonRP(radius, maximumRadius, diffuseness);
00080 } else if(A <= 19 && A > 6) {
00081 const G4double radius = ParticleTable::getRadiusParameter(A, Z);
00082 const G4double diffuseness = ParticleTable::getSurfaceDiffuseness(A, Z);
00083 const G4double maximumRadius = ParticleTable::getMaximumNuclearRadius(A, Z);
00084 rpCorrelationFunction = new NuclearDensityFunctions::ModifiedHarmonicOscillatorRP(radius, maximumRadius, diffuseness);
00085 } else if(A <= 6 && A > 1) {
00086 const G4double radius = ParticleTable::getRadiusParameter(A, Z);
00087 const G4double maximumRadius = ParticleTable::getMaximumNuclearRadius(A, Z);
00088 rpCorrelationFunction = new NuclearDensityFunctions::GaussianRP(maximumRadius, Math::oneOverSqrtThree * radius);
00089 } else {
00090 ERROR("No r-p correlation function for target A = "
00091 << A << " Z = " << Z << std::endl);
00092 return NULL;
00093 }
00094
00095 class InverseCDFOneThird : public IFunction1D {
00096 public:
00097 InverseCDFOneThird(IFunction1D const * const f) :
00098 IFunction1D(f->getXMinimum(), f->getXMaximum()),
00099 theFunction(f),
00100 normalisation(1./theFunction->integrate(xMin,xMax))
00101 {}
00102
00103 G4double operator()(const G4double x) const {
00104 return Math::pow13(normalisation * theFunction->integrate(xMin,x));
00105 }
00106 private:
00107 IFunction1D const * const theFunction;
00108 const G4double normalisation;
00109 } *theInverseCDFOneThird = new InverseCDFOneThird(rpCorrelationFunction);
00110
00111 InverseInterpolationTable *theTable = new InverseInterpolationTable(*theInverseCDFOneThird);
00112 delete theInverseCDFOneThird;
00113 delete rpCorrelationFunction;
00114 DEBUG("Creating r-p correlation function for A=" << A << ", Z=" << Z << ":"
00115 << std::endl << theTable->print() << std::endl);
00116
00117 rpCorrelationTableCache[nuclideID] = theTable;
00118 return theTable;
00119 } else {
00120 return mapEntry->second;
00121 }
00122 }
00123
00124 InverseInterpolationTable *NuclearDensityFactory::createRCDFTable(const G4int A, const G4int Z) {
00125 const G4int nuclideID = 1000*Z + A;
00126 const std::map<G4int,InverseInterpolationTable*>::const_iterator mapEntry = rCDFTableCache.find(nuclideID);
00127 if(mapEntry == rCDFTableCache.end()) {
00128
00129 IFunction1D *rDensityFunction;
00130 if(A > 19) {
00131 G4double radius = ParticleTable::getRadiusParameter(A, Z);
00132 G4double diffuseness = ParticleTable::getSurfaceDiffuseness(A, Z);
00133 G4double maximumRadius = ParticleTable::getMaximumNuclearRadius(A, Z);
00134 rDensityFunction = new NuclearDensityFunctions::WoodsSaxon(radius, maximumRadius, diffuseness);
00135 } else if(A <= 19 && A > 6) {
00136 G4double radius = ParticleTable::getRadiusParameter(A, Z);
00137 G4double diffuseness = ParticleTable::getSurfaceDiffuseness(A, Z);
00138 G4double maximumRadius = ParticleTable::getMaximumNuclearRadius(A, Z);
00139 rDensityFunction = new NuclearDensityFunctions::ModifiedHarmonicOscillator(radius, maximumRadius, diffuseness);
00140 } else if(A <= 6 && A > 2) {
00141 G4double radius = ParticleTable::getRadiusParameter(A, Z);
00142 G4double maximumRadius = ParticleTable::getMaximumNuclearRadius(A, Z);
00143 rDensityFunction = new NuclearDensityFunctions::Gaussian(maximumRadius, Math::oneOverSqrtThree * radius);
00144 } else if(A == 2 && Z == 1) {
00145 rDensityFunction = new NuclearDensityFunctions::ParisR();
00146 } else {
00147 ERROR("No nuclear density function for target A = "
00148 << A << " Z = " << Z << std::endl);
00149 return NULL;
00150 }
00151
00152 InverseInterpolationTable *theTable = rDensityFunction->inverseCDFTable();
00153 delete rDensityFunction;
00154 DEBUG("Creating inverse position CDF for A=" << A << ", Z=" << Z << ":" <<
00155 std::endl << theTable->print() << std::endl);
00156
00157 rCDFTableCache[nuclideID] = theTable;
00158 return theTable;
00159 } else {
00160 return mapEntry->second;
00161 }
00162 }
00163
00164 InverseInterpolationTable *NuclearDensityFactory::createPCDFTable(const G4int A, const G4int Z) {
00165 const G4int nuclideID = 1000*Z + A;
00166 const std::map<G4int,InverseInterpolationTable*>::const_iterator mapEntry = pCDFTableCache.find(nuclideID);
00167 if(mapEntry == pCDFTableCache.end()) {
00168 IFunction1D *pDensityFunction;
00169 if(A > 19) {
00170 pDensityFunction = new NuclearDensityFunctions::HardSphere(PhysicalConstants::Pf);
00171 } else if(A <= 19 && A > 2) {
00172 G4double momentumRMS = Math::oneOverSqrtThree * ParticleTable::getMomentumRMS(A, Z);
00173 pDensityFunction = new NuclearDensityFunctions::Gaussian(5.*momentumRMS, momentumRMS);
00174 } else if(A == 2 && Z == 1) {
00175 pDensityFunction = new NuclearDensityFunctions::ParisP();
00176 } else {
00177 ERROR("No nuclear density function for target A = "
00178 << A << " Z = " << Z << std::endl);
00179 return NULL;
00180 }
00181
00182 InverseInterpolationTable *theTable = pDensityFunction->inverseCDFTable();
00183 delete pDensityFunction;
00184 DEBUG("Creating inverse momentum CDF for A=" << A << ", Z=" << Z << ":" <<
00185 std::endl << theTable->print() << std::endl);
00186
00187 pCDFTableCache[nuclideID] = theTable;
00188 return theTable;
00189 } else {
00190 return mapEntry->second;
00191 }
00192 }
00193
00194 ParticleSampler *NuclearDensityFactory::createParticleSampler(const G4int A, const G4int Z) {
00195 InverseInterpolationTable *rCDFTable = NuclearDensityFactory::createRCDFTable(A, Z);
00196 InverseInterpolationTable *pCDFTable = NuclearDensityFactory::createPCDFTable(A, Z);
00197 return new ParticleSampler(A, Z, rCDFTable, pCDFTable);
00198 }
00199
00200 }