ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4INCLParticleSampler.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4INCLParticleSampler.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 
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));
64  }
65 
67  }
68 
70  theDensity = d;
72  }
73 
75  thePotential = p;
77  }
78 
80  if(theDensity && thePotential) {
81  if(rpCorrelationCoefficient[Proton]>0.99999) {
83  } else {
85  }
86  if(rpCorrelationCoefficient[Neutron]>0.99999) {
88  } else {
90  }
91  } else {
94  }
95  }
96 
98  ParticleList aList;
99  sampleParticlesIntoList(position, aList);
100  return aList;
101  }
102 
104 
106  // sampling without correlation, we need to initialize the CDF tables
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 
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 
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 
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