Geant4_10
G4INCLSurfaceAvatar.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 
37 /*
38  * G4INCLReflectionAvatar.cc
39  *
40  * \date Jun 8, 2009
41  * \author Pekka Kaitaniemi
42  */
43 
44 #include "G4INCLSurfaceAvatar.hh"
45 #include "G4INCLRandom.hh"
48 #include "G4INCLClustering.hh"
49 #include <sstream>
50 #include <string>
51 
52 namespace G4INCL {
53 
55  :IAvatar(time), theParticle(aParticle), theNucleus(n),
56  particlePIn(0.),
57  particlePOut(0.),
58  particleTOut(0.),
59  TMinusV(0.),
60  TMinusV2(0.),
61  particleMass(0.),
62  sinIncidentAngle(0.),
63  cosIncidentAngle(0.),
64  sinRefractionAngle(0.),
65  cosRefractionAngle(0.),
66  refractionIndexRatio(0.),
67  internalReflection(false)
68  {
70  }
71 
73 
74  }
75 
77  if(theParticle->isTargetSpectator()) {
78  INCL_DEBUG("Particle " << theParticle->getID() << " is a spectator, reflection" << std::endl);
79  return new ReflectionChannel(theNucleus, theParticle);
80  }
81 
82  // We forbid transmission of resonances below the Fermi energy. Emitting a
83  // delta particle below Tf can lead to negative excitation energies, since
84  // CDPP assumes that particles stay in the Fermi sea.
85  if(theParticle->isResonance()) {
86  const G4double theFermiEnergy = theNucleus->getPotential()->getFermiEnergy(theParticle);
87  if(theParticle->getKineticEnergy()<theFermiEnergy) {
88  INCL_DEBUG("Particle " << theParticle->getID() << " is a resonance below Tf, reflection" << std::endl
89  << " Tf=" << theFermiEnergy << ", EKin=" << theParticle->getKineticEnergy() << std::endl);
90  return new ReflectionChannel(theNucleus, theParticle);
91  }
92  }
93 
94  // Don't try to make a cluster if the leading particle is too slow
95  const G4double transmissionProbability = getTransmissionProbability(theParticle);
96  const G4double TOut = TMinusV;
97  const G4double kOut = particlePOut;
98  const G4double cosR = cosRefractionAngle;
99 
100  INCL_DEBUG("Transmission probability for particle " << theParticle->getID() << " = " << transmissionProbability << std::endl);
101  /* Don't attempt to construct clusters when a projectile spectator is
102  * trying to escape during a nucleus-nucleus collision. The idea behind
103  * this is that projectile spectators will later be collected in the
104  * projectile remnant, and trying to clusterise them somewhat feels like
105  * G4double counting. Moreover, applying the clustering algorithm on escaping
106  * projectile spectators makes the code *really* slow if the projectile is
107  * large.
108  */
109  if(theParticle->isNucleon()
110  && (!theParticle->isProjectileSpectator() || !theNucleus->isNucleusNucleusCollision())
111  && transmissionProbability>1.E-4) {
112  Cluster *candidateCluster = 0;
113 
114  candidateCluster = Clustering::getCluster(theNucleus, theParticle);
115  if(candidateCluster != 0 &&
116  Clustering::clusterCanEscape(theNucleus, candidateCluster)) {
117 
118  INCL_DEBUG("Cluster algorithm succeded. Candidate cluster:" << std::endl << candidateCluster->print() << std::endl);
119 
120  // Check if the cluster can penetrate the Coulomb barrier
121  const G4double clusterTransmissionProbability = getTransmissionProbability(candidateCluster);
122  const G4double x = Random::shoot();
123 
124  INCL_DEBUG("Transmission probability for cluster " << candidateCluster->getID() << " = " << clusterTransmissionProbability << std::endl);
125 
126  if (x <= clusterTransmissionProbability) {
127  theNucleus->getStore()->getBook().incrementEmittedClusters();
128  INCL_DEBUG("Cluster " << candidateCluster->getID() << " passes the Coulomb barrier, transmitting." << std::endl);
129  return new TransmissionChannel(theNucleus, candidateCluster);
130  } else {
131  INCL_DEBUG("Cluster " << candidateCluster->getID() << " does not pass the Coulomb barrier. Falling back to transmission of the leading particle." << std::endl);
132  delete candidateCluster;
133  }
134  } else {
135  delete candidateCluster;
136  }
137  }
138 
139  // If we haven't transmitted a cluster (maybe cluster feature was
140  // disabled or maybe we just can't produce an acceptable cluster):
141 
142  // Always transmit projectile spectators if no cluster was formed and if
143  // transmission is energetically allowed
144  if(theParticle->isProjectileSpectator() && transmissionProbability>0.) {
145  INCL_DEBUG("Particle " << theParticle->getID() << " is a projectile spectator, transmission" << std::endl);
146  return new TransmissionChannel(theNucleus, theParticle, TOut);
147  }
148 
149  // Transmit or reflect depending on the transmission probability
150  const G4double x = Random::shoot();
151 
152  if(x <= transmissionProbability) { // Transmission
153  INCL_DEBUG("Particle " << theParticle->getID() << " passes the Coulomb barrier, transmitting." << std::endl);
154  if(theNucleus->getStore()->getConfig()->getRefraction()) {
155  return new TransmissionChannel(theNucleus, theParticle, kOut, cosR);
156  } else {
157  return new TransmissionChannel(theNucleus, theParticle, TOut);
158  }
159  } else { // Reflection
160  INCL_DEBUG("Particle " << theParticle->getID() << " does not pass the Coulomb barrier, reflection." << std::endl);
161  return new ReflectionChannel(theNucleus, theParticle);
162  }
163  }
164 
166  return getChannel()->getFinalState();
167  }
168 
170 
172  ParticleList outgoing = fs->getOutgoingParticles();
173  if(!outgoing.empty()) { // Transmission
174 // assert(outgoing.size()==1);
175  Particle *out = outgoing.front();
176  out->rpCorrelate();
177  if(out->isCluster()) {
178  Cluster *clusterOut = dynamic_cast<Cluster*>(out);
179  ParticleList const components = clusterOut->getParticles();
180  for(ParticleIter i=components.begin(), e=components.end(); i!=e; ++i) {
181  if(!(*i)->isTargetSpectator())
182  theNucleus->getStore()->getBook().decrementCascading();
183  }
184  } else if(!theParticle->isTargetSpectator()) {
185 // assert(out==theParticle);
186  theNucleus->getStore()->getBook().decrementCascading();
187  }
188  }
189  return fs;
190  }
191 
192  std::string SurfaceAvatar::dump() const {
193  std::stringstream ss;
194  ss << "(avatar " << theTime << " 'reflection" << std::endl
195  << "(list " << std::endl
196  << theParticle->dump()
197  << "))" << std::endl;
198  return ss.str();
199  }
200 
202 
203  particleMass = particle->getMass();
204  const G4double V = particle->getPotentialEnergy();
205 
206  // Correction to the particle kinetic energy if using real masses
207  const G4int theA = theNucleus->getA();
208  const G4int theZ = theNucleus->getZ();
209  const G4double correction = particle->getEmissionQValueCorrection(theA, theZ);
210  particleTOut = particle->getKineticEnergy() + correction;
211 
212  if (particleTOut <= V) // No transmission if total energy < 0
213  return 0.0;
214 
215  TMinusV = particleTOut-V;
216  TMinusV2 = TMinusV*TMinusV;
217 
218  // Momenta in and out
219  const G4double particlePIn2 = particle->getMomentum().mag2();
220  const G4double particlePOut2 = 2.*particleMass*TMinusV+TMinusV2;
221  particlePIn = std::sqrt(particlePIn2);
222  particlePOut = std::sqrt(particlePOut2);
223 
224  // Compute the transmission probability
225  G4double theTransmissionProbability;
226  if(theNucleus->getStore()->getConfig()->getRefraction()) {
227  // Use the formula with refraction
228  initializeRefractionVariables(particle);
229 
230  if(internalReflection)
231  return 0.; // total internal reflection
232 
233  // Intermediate variables for calculation
234  const G4double x = refractionIndexRatio*cosIncidentAngle;
235  const G4double y = (x - cosRefractionAngle) / (x + cosRefractionAngle);
236 
237  theTransmissionProbability = 1. - y*y;
238  } else {
239  // Use the formula without refraction
240 
241  // Intermediate variable for calculation
242  const G4double y = particlePIn+particlePOut;
243 
244  // The transmission probability for a potential step
245  theTransmissionProbability = 4.*particlePIn*particlePOut/(y*y);
246  }
247 
248  // For neutral and negative particles, no Coulomb transmission
249  // Also, no Coulomb if the particle takes away all of the nuclear charge
250  const G4int particleZ = particle->getZ();
251  if (particleZ <= 0 || particleZ >= theZ)
252  return theTransmissionProbability;
253 
254  // Nominal Coulomb barrier
255  const G4double theTransmissionBarrier = theNucleus->getTransmissionBarrier(particle);
256  if (TMinusV >= theTransmissionBarrier) // Above the Coulomb barrier
257  return theTransmissionProbability;
258 
259  // Coulomb-penetration factor
260  const G4double px = std::sqrt(TMinusV/theTransmissionBarrier);
261  const G4double logCoulombTransmission =
262  particleZ*(theZ-particleZ)/137.03*std::sqrt(2.*particleMass/TMinusV/(1.+TMinusV/2./particleMass))
263  *(std::acos(px)-px*std::sqrt(1.-px*px));
264  INCL_DEBUG("Coulomb barrier, logCoulombTransmission=" << logCoulombTransmission << std::endl);
265  if (logCoulombTransmission > 35.) // Transmission is forbidden by Coulomb
266  return 0.;
267  theTransmissionProbability *= std::exp(-2.*logCoulombTransmission);
268 
269  return theTransmissionProbability;
270  }
271 
272  void SurfaceAvatar::initializeRefractionVariables(Particle const * const particle) {
273  cosIncidentAngle = particle->getCosRPAngle();
274  if(cosIncidentAngle>1.)
275  cosIncidentAngle=1.;
276  sinIncidentAngle = std::sqrt(1. - cosIncidentAngle*cosIncidentAngle);
277  refractionIndexRatio = particlePIn/particlePOut;
278  const G4double sinCandidate = refractionIndexRatio*sinIncidentAngle;
279  internalReflection = (std::fabs(sinCandidate)>1.);
280  if(internalReflection) {
281  sinRefractionAngle = 1.;
282  cosRefractionAngle = 0.;
283  } else {
284  sinRefractionAngle = sinCandidate;
285  cosRefractionAngle = std::sqrt(1. - sinRefractionAngle*sinRefractionAngle);
286  }
287  INCL_DEBUG("Refraction parameters initialised as follows:\n"
288  << " cosIncidentAngle=" << cosIncidentAngle << '\n'
289  << " sinIncidentAngle=" << sinIncidentAngle << '\n'
290  << " cosRefractionAngle=" << cosRefractionAngle << '\n'
291  << " sinRefractionAngle=" << sinRefractionAngle << '\n'
292  << " refractionIndexRatio=" << refractionIndexRatio << '\n'
293  << " internalReflection=" << internalReflection << std::endl);
294  }
295 }
G4int getA() const
Returns the baryon number.
G4double getCosRPAngle() const
Get the cosine of the angle between position and momentum.
G4bool isResonance() const
Is it a resonance?
G4INCL::IChannel * getChannel()
G4double getMass() const
Get the cached particle mass.
virtual FinalState * postInteraction(FinalState *)
ParticleList const & getParticles() const
G4bool isNucleusNucleusCollision() const
Is it a nucleus-nucleus collision?
G4INCL::FinalState * getFinalState()
G4bool isTargetSpectator() const
const G4INCL::ThreeVector & getMomentum() const
Config const * getConfig()
Definition: G4INCLStore.hh:251
Store * getStore() const
tuple x
Definition: test.py:50
G4double getTransmissionProbability(Particle const *const particle)
Calculate the transmission probability for the particle.
int G4int
Definition: G4Types.hh:78
G4double mag2() const
G4double getPotentialEnergy() const
Get the particle potential energy.
static G4bool clusterCanEscape(Nucleus const *const n, Cluster const *const c)
G4double getFermiEnergy(const Particle *const p) const
Return the Fermi energy for a particle.
void rpCorrelate()
Make the particle follow a strict r-p correlation.
ParticleList const & getOutgoingParticles() const
Double_t y
Definition: plot.C:279
Char_t n[5]
void setType(AvatarType t)
SurfaceAvatar(G4INCL::Particle *aParticle, G4double time, G4INCL::Nucleus *aNucleus)
Book & getBook()
Definition: G4INCLStore.hh:237
std::string print() const
G4int getZ() const
Returns the charge number.
void incrementEmittedClusters()
Definition: G4INCLBook.hh:78
NuclearPotential::INuclearPotential const * getPotential() const
Getter for thePotential.
void decrementCascading()
Definition: G4INCLBook.hh:77
Static class for cluster formation.
G4double getTransmissionBarrier(Particle const *const p)
Get the transmission barrier.
static Cluster * getCluster(Nucleus *n, Particle *p)
G4bool isProjectileSpectator() const
virtual G4INCL::FinalState * getFinalState()=0
std::string dump() const
G4bool isNucleon() const
G4bool isCluster() const
G4double shoot()
Definition: G4INCLRandom.cc:74
G4double getKineticEnergy() const
Get the particle kinetic energy.
double G4double
Definition: G4Types.hh:76
G4double getEmissionQValueCorrection(const G4int AParent, const G4int ZParent) const
Computes correction on the emission Q-value.
#define INCL_DEBUG(x)
std::string dump() const
ParticleList::const_iterator ParticleIter
long getID() const
G4bool getRefraction() const
True if we should use refraction.