Geant4_10
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 // Pekka Kaitaniemi, CEA and Helsinki Institute of Physics
28 // Davide Mancusi, CEA
29 // Alain Boudard, CEA
30 // Sylvie Leray, CEA
31 // Joseph Cugnon, University of Liege
32 //
33 #define INCLXX_IN_GEANT4_MODE 1
34 
35 #include "globals.hh"
36 
44 #include "G4INCLParticleSampler.hh"
45 #include "G4INCLParticleTable.hh"
47 
48 namespace G4INCL {
49 
51  sampleOneProton(&ParticleSampler::sampleOneParticleWithoutRPCorrelation),
52  sampleOneNeutron(&ParticleSampler::sampleOneParticleWithoutRPCorrelation),
53  theA(A),
54  theZ(Z),
55  theDensity(NULL),
56  thePotential(NULL)
57  {
58  std::fill(theRCDFTable, theRCDFTable + UnknownParticle, static_cast<InverseInterpolationTable *>(NULL));
59  std::fill(thePCDFTable, thePCDFTable + UnknownParticle, static_cast<InverseInterpolationTable *>(NULL));
60  std::fill(rpCorrelationCoefficient, rpCorrelationCoefficient + UnknownParticle, 1.);
61  rpCorrelationCoefficient[Proton] = ParticleTable::getRPCorrelationCoefficient(Proton);
62  rpCorrelationCoefficient[Neutron] = ParticleTable::getRPCorrelationCoefficient(Neutron);
63  }
64 
66  }
67 
69  theDensity = d;
70  updateSampleOneParticleMethods();
71  }
72 
74  thePotential = p;
75  updateSampleOneParticleMethods();
76  }
77 
78  void ParticleSampler::updateSampleOneParticleMethods() {
79  if(theDensity && thePotential) {
80  if(rpCorrelationCoefficient[Proton]>0.99999) {
81  sampleOneProton = &ParticleSampler::sampleOneParticleWithRPCorrelation;
82  } else {
83  sampleOneProton = &ParticleSampler::sampleOneParticleWithFuzzyRPCorrelation;
84  }
85  if(rpCorrelationCoefficient[Neutron]>0.99999) {
86  sampleOneNeutron = &ParticleSampler::sampleOneParticleWithRPCorrelation;
87  } else {
88  sampleOneNeutron = &ParticleSampler::sampleOneParticleWithFuzzyRPCorrelation;
89  }
90  } else {
91  sampleOneProton = &ParticleSampler::sampleOneParticleWithoutRPCorrelation;
92  sampleOneNeutron = &ParticleSampler::sampleOneParticleWithoutRPCorrelation;
93  }
94  }
95 
97 
98  if(sampleOneProton == &ParticleSampler::sampleOneParticleWithoutRPCorrelation) {
99  // sampling without correlation, we need to initialize the CDF tables
100  theRCDFTable[Proton] = NuclearDensityFactory::createRCDFTable(Proton, theA, theZ);
101  thePCDFTable[Proton] = NuclearDensityFactory::createPCDFTable(Proton, theA, theZ);
102  theRCDFTable[Neutron] = NuclearDensityFactory::createRCDFTable(Neutron, theA, theZ);
103  thePCDFTable[Neutron] = NuclearDensityFactory::createPCDFTable(Neutron, theA, theZ);
104  }
105 
106  ParticleList theList;
107  if(theA > 2) {
108  ParticleType type = Proton;
109  ParticleSamplerMethod sampleOneParticle = sampleOneProton;
110  for(G4int i = 1; i <= theA; ++i) {
111  if(i == (theZ + 1)) { // Nucleons [Z+1..A] are neutrons
112  type = Neutron;
113  sampleOneParticle = sampleOneNeutron;
114  }
115  Particle *p = (this->*sampleOneParticle)(type);
116  p->setPosition(position + p->getPosition());
117  theList.push_back(p);
118  }
119  } else {
120  // For deuterons, only sample the proton position and momentum. The
121  // neutron position and momenta are determined by the conditions of
122  // vanishing CM position and total momentum.
123 // assert(theZ==1);
124  Particle *aProton = (this->*(this->sampleOneProton))(Proton);
125  Particle *aNeutron = new Particle(Neutron, -aProton->getMomentum(), position - aProton->getPosition());
126  aProton->setPosition(position + aProton->getPosition());
127  theList.push_back(aProton);
128  theList.push_back(aNeutron);
129  }
130 
131  return theList;
132  }
133 
134  Particle *ParticleSampler::sampleOneParticleWithRPCorrelation(const ParticleType t) const {
135 // assert(theDensity && thePotential);
136  const G4double theFermiMomentum = thePotential->getFermiMomentum(t);
137  const ThreeVector momentumVector = Random::sphereVector(theFermiMomentum);
138  const G4double momentumAbs = momentumVector.mag();
139  const G4double momentumRatio = momentumAbs/theFermiMomentum;
140  const G4double reflectionRadius = theDensity->getMaxRFromP(t, momentumRatio);
141  const ThreeVector positionVector = Random::sphereVector(reflectionRadius);
142  Particle *aParticle = new Particle(t, momentumVector, positionVector);
143  aParticle->setUncorrelatedMomentum(momentumAbs);
144  return aParticle;
145  }
146 
147  Particle *ParticleSampler::sampleOneParticleWithoutRPCorrelation(const ParticleType t) const {
148  const G4double position = (*(theRCDFTable[t]))(Random::shoot());
149  const G4double momentum = (*(thePCDFTable[t]))(Random::shoot());
150  ThreeVector positionVector = Random::normVector(position);
151  ThreeVector momentumVector = Random::normVector(momentum);
152  return new Particle(t, momentumVector, positionVector);
153  }
154 
155  Particle *ParticleSampler::sampleOneParticleWithFuzzyRPCorrelation(const ParticleType t) const {
156 // assert(theDensity && thePotential);
157  std::pair<G4double,G4double> ranNumbers = Random::correlatedUniform(rpCorrelationCoefficient[t]);
158  const G4double x = Math::pow13(ranNumbers.first);
159  const G4double y = Math::pow13(ranNumbers.second);
160  const G4double theFermiMomentum = thePotential->getFermiMomentum(t);
161  const ThreeVector momentumVector = Random::normVector(y*theFermiMomentum);
162  const G4double reflectionRadius = theDensity->getMaxRFromP(t, x);
163  const ThreeVector positionVector = Random::sphereVector(reflectionRadius);
164  Particle *aParticle = new Particle(t, momentumVector, positionVector);
165  aParticle->setUncorrelatedMomentum(x*theFermiMomentum);
166  return aParticle;
167  }
168 
169 }
170 
ParticleList sampleParticles(ThreeVector const &position)
ThreeVector sphereVector(G4double rmax=1.)
Float_t d
Definition: plot.C:237
G4double getFermiMomentum(const Particle *const p) const
Return the Fermi momentum for a particle.
InverseInterpolationTable * createPCDFTable(const ParticleType t, const G4int A, const G4int Z)
const char * p
Definition: xmltok.h:285
const G4INCL::ThreeVector & getMomentum() const
tuple x
Definition: test.py:50
void setUncorrelatedMomentum(const G4double p)
Set the uncorrelated momentum.
int G4int
Definition: G4Types.hh:78
ThreeVector normVector(G4double norm=1.)
Double_t y
Definition: plot.C:279
Class for sampling particles in a nucleus.
void setDensity(NuclearDensity const *const d)
Setter for theDensity.
Float_t Z
Definition: plot.C:39
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.
InverseInterpolationTable * createRCDFTable(const ParticleType t, const G4int A, const G4int Z)
G4double shoot()
Definition: G4INCLRandom.cc:74
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.