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
00044 #include "G4INCLParticleSampler.hh"
00045
00046 namespace G4INCL {
00047
00048 ParticleSampler::ParticleSampler(const G4int A, const G4int Z, InverseInterpolationTable const * const rCDFTable, InverseInterpolationTable const * const pCDFTable) :
00049 sampleOneParticle(&ParticleSampler::sampleOneParticleWithoutRPCorrelation),
00050 theA(A),
00051 theZ(Z),
00052 theRCDFTable(rCDFTable),
00053 thePCDFTable(pCDFTable),
00054 theDensity(NULL),
00055 thePotential(NULL)
00056 {}
00057
00058 ParticleSampler::~ParticleSampler() {
00059 }
00060
00061 void ParticleSampler::setDensity(NuclearDensity const * const d) {
00062 theDensity = d;
00063 updateSampleOneParticleMethod();
00064 }
00065
00066 void ParticleSampler::setPotential(NuclearPotential::INuclearPotential const * const p) {
00067 thePotential = p;
00068 updateSampleOneParticleMethod();
00069 }
00070
00071 void ParticleSampler::updateSampleOneParticleMethod() {
00072 if(theDensity && thePotential)
00073 sampleOneParticle = &ParticleSampler::sampleOneParticleWithRPCorrelation;
00074 else
00075 sampleOneParticle = &ParticleSampler::sampleOneParticleWithoutRPCorrelation;
00076 }
00077
00078 ParticleList ParticleSampler::sampleParticles(const ThreeVector &position) const {
00079 ParticleList theList;
00080 if(theA > 2) {
00081 ParticleType type = Proton;
00082 for(G4int i = 1; i <= theA; ++i) {
00083 if(i == (theZ + 1))
00084 type = Neutron;
00085 Particle *p = (this->*(this->sampleOneParticle))(type);
00086 p->setPosition(position + p->getPosition());
00087 theList.push_back(p);
00088 }
00089 } else {
00090
00091
00092
00093
00094 Particle *aProton = (this->*(this->sampleOneParticle))(Proton);
00095 Particle *aNeutron = new Particle(Neutron, -aProton->getMomentum(), position - aProton->getPosition());
00096 aProton->setPosition(position + aProton->getPosition());
00097 theList.push_back(aProton);
00098 theList.push_back(aNeutron);
00099 }
00100
00101 return theList;
00102 }
00103
00104 Particle *ParticleSampler::sampleOneParticleWithRPCorrelation(const ParticleType t) const {
00105
00106 const G4double theFermiMomentum = thePotential->getFermiMomentum(t);
00107 const ThreeVector momentumVector = Random::sphereVector(theFermiMomentum);
00108 const G4double momentumRatio = momentumVector.mag()/theFermiMomentum;
00109 const ThreeVector positionVector = Random::sphereVector(theDensity->getMaxRFromP(momentumRatio));
00110 return new Particle(t, momentumVector, positionVector);
00111 }
00112
00113 Particle *ParticleSampler::sampleOneParticleWithoutRPCorrelation(const ParticleType t) const {
00114 const G4double position = (*theRCDFTable)(Random::shoot());
00115 const G4double momentum = (*thePCDFTable)(Random::shoot());
00116 ThreeVector positionVector = Random::normVector(position);
00117 ThreeVector momentumVector = Random::normVector(momentum);
00118 return new Particle(t, momentumVector, positionVector);
00119 }
00120
00121 }
00122