Geant4  10.03.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4INCLParticleSampler.cc
Go to the documentation of this file.
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 
45 #include "G4INCLParticleSampler.hh"
46 #include "G4INCLParticleTable.hh"
48 
49 namespace G4INCL {
50 
52  sampleOneProton(&ParticleSampler::sampleOneParticleWithoutRPCorrelation),
53  sampleOneNeutron(&ParticleSampler::sampleOneParticleWithoutRPCorrelation),
54  theA(A),
55  theZ(Z),
56  theDensity(NULL),
57  thePotential(NULL)
58  {
59  std::fill(theRCDFTable, theRCDFTable + UnknownParticle, static_cast<InterpolationTable *>(NULL));
60  std::fill(thePCDFTable, thePCDFTable + UnknownParticle, static_cast<InterpolationTable *>(NULL));
61  std::fill(rpCorrelationCoefficient, rpCorrelationCoefficient + UnknownParticle, 1.);
62  rpCorrelationCoefficient[Proton] = ParticleTable::getRPCorrelationCoefficient(Proton);
63  rpCorrelationCoefficient[Neutron] = ParticleTable::getRPCorrelationCoefficient(Neutron);
64  }
65 
67  }
68 
70  theDensity = d;
71  updateSampleOneParticleMethods();
72  }
73 
75  thePotential = p;
76  updateSampleOneParticleMethods();
77  }
78 
79  void ParticleSampler::updateSampleOneParticleMethods() {
80  if(theDensity && thePotential) {
81  if(rpCorrelationCoefficient[Proton]>0.99999) {
82  sampleOneProton = &ParticleSampler::sampleOneParticleWithRPCorrelation;
83  } else {
84  sampleOneProton = &ParticleSampler::sampleOneParticleWithFuzzyRPCorrelation;
85  }
86  if(rpCorrelationCoefficient[Neutron]>0.99999) {
87  sampleOneNeutron = &ParticleSampler::sampleOneParticleWithRPCorrelation;
88  } else {
89  sampleOneNeutron = &ParticleSampler::sampleOneParticleWithFuzzyRPCorrelation;
90  }
91  } else {
92  sampleOneProton = &ParticleSampler::sampleOneParticleWithoutRPCorrelation;
93  sampleOneNeutron = &ParticleSampler::sampleOneParticleWithoutRPCorrelation;
94  }
95  }
96 
98  ParticleList aList;
99  sampleParticlesIntoList(position, aList);
100  return aList;
101  }
102 
104 
105  if(sampleOneProton == &ParticleSampler::sampleOneParticleWithoutRPCorrelation) {
106  // sampling without correlation, we need to initialize the CDF tables
107  theRCDFTable[Proton] = NuclearDensityFactory::createRCDFTable(Proton, theA, theZ);
108  thePCDFTable[Proton] = NuclearDensityFactory::createPCDFTable(Proton, theA, theZ);
109  theRCDFTable[Neutron] = NuclearDensityFactory::createRCDFTable(Neutron, theA, theZ);
110  thePCDFTable[Neutron] = NuclearDensityFactory::createPCDFTable(Neutron, theA, theZ);
111  }
112 
113  theList.resize(theA);
114  if(theA > 2) {
115  ParticleType type = Proton;
116  ParticleSamplerMethod sampleOneParticle = sampleOneProton;
117  for(G4int i = 0; i < theA; ++i) {
118  if(i == theZ) { // Nucleons [Z..A-1] are neutrons
119  type = Neutron;
120  sampleOneParticle = sampleOneNeutron;
121  }
122  Particle *p = (this->*sampleOneParticle)(type);
123  p->setPosition(position + p->getPosition());
124  theList[i] = p;
125  }
126  } else {
127  // For deuterons, only sample the proton position and momentum. The
128  // neutron position and momenta are determined by the conditions of
129  // vanishing CM position and total momentum.
130 // assert(theZ==1);
131  Particle *aProton = (this->*(this->sampleOneProton))(Proton);
132  Particle *aNeutron = new Particle(Neutron, -aProton->getMomentum(), position - aProton->getPosition());
133  aProton->setPosition(position + aProton->getPosition());
134  theList[0] = aProton;
135  theList[1] = aNeutron;
136  }
137  }
138 
139  Particle *ParticleSampler::sampleOneParticleWithRPCorrelation(const ParticleType t) const {
140 // assert(theDensity && thePotential);
141  const G4double theFermiMomentum = thePotential->getFermiMomentum(t);
142  const ThreeVector momentumVector = Random::sphereVector(theFermiMomentum);
143  const G4double momentumAbs = momentumVector.mag();
144  const G4double momentumRatio = momentumAbs/theFermiMomentum;
145  const G4double reflectionRadius = theDensity->getMaxRFromP(t, momentumRatio);
146  const ThreeVector positionVector = Random::sphereVector(reflectionRadius);
147  Particle *aParticle = new Particle(t, momentumVector, positionVector);
148  aParticle->setUncorrelatedMomentum(momentumAbs);
149  return aParticle;
150  }
151 
152  Particle *ParticleSampler::sampleOneParticleWithoutRPCorrelation(const ParticleType t) const {
153  const G4double position = (*(theRCDFTable[t]))(Random::shoot());
154  const G4double momentum = (*(thePCDFTable[t]))(Random::shoot());
155  ThreeVector positionVector = Random::normVector(position);
156  ThreeVector momentumVector = Random::normVector(momentum);
157  return new Particle(t, momentumVector, positionVector);
158  }
159 
160  Particle *ParticleSampler::sampleOneParticleWithFuzzyRPCorrelation(const ParticleType t) const {
161 // assert(theDensity && thePotential);
162  std::pair<G4double,G4double> ranNumbers = Random::correlatedUniform(rpCorrelationCoefficient[t]);
163  const G4double x = Math::pow13(ranNumbers.first);
164  const G4double y = Math::pow13(ranNumbers.second);
165  const G4double theFermiMomentum = thePotential->getFermiMomentum(t);
166  const ThreeVector momentumVector = Random::normVector(y*theFermiMomentum);
167  const G4double reflectionRadius = theDensity->getMaxRFromP(t, x);
168  const ThreeVector positionVector = Random::sphereVector(reflectionRadius);
169  Particle *aParticle = new Particle(t, momentumVector, positionVector);
170  aParticle->setUncorrelatedMomentum(x*theFermiMomentum);
171  return aParticle;
172  }
173 
174 }
175 
ParticleList sampleParticles(ThreeVector const &position)
ThreeVector sphereVector(G4double rmax=1.)
InterpolationTable * createPCDFTable(const ParticleType t, const G4int A, const G4int Z)
G4double getFermiMomentum(const Particle *const p) const
Return the Fermi momentum for a particle.
const char * p
Definition: xmltok.h:285
const G4INCL::ThreeVector & getMomentum() const
tuple x
Definition: test.py:50
InterpolationTable * createRCDFTable(const ParticleType t, const G4int A, const G4int Z)
void setUncorrelatedMomentum(const G4double p)
Set the uncorrelated momentum.
int G4int
Definition: G4Types.hh:78
ThreeVector normVector(G4double norm=1.)
void sampleParticlesIntoList(ThreeVector const &position, ParticleList &theList)
Class for sampling particles in a nucleus.
double A(double temperature)
void setDensity(NuclearDensity const *const d)
Setter for theDensity.
std::pair< G4double, G4double > correlatedUniform(const G4double corrCoeff)
Generate pairs of correlated uniform random numbers.
virtual void setPosition(const G4INCL::ThreeVector &position)
ParticleSampler(const G4int A, const G4int Z)
Constructor.
const G4INCL::ThreeVector & getPosition() const
void setPotential(NuclearPotential::INuclearPotential const *const p)
Setter for thePotential.
G4double shoot()
Definition: G4INCLRandom.cc:93
G4double mag() const
G4double getRPCorrelationCoefficient(const ParticleType t)
Get the value of the r-p correlation coefficient.
double G4double
Definition: G4Types.hh:76
G4double pow13(G4double x)
G4double getMaxRFromP(const ParticleType t, const G4double p) const
Get the maximum allowed radius for a given momentum.