Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 // 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 
38 /*
39  * G4INCLReflectionAvatar.cc
40  *
41  * \date Jun 8, 2009
42  * \author Pekka Kaitaniemi
43  */
44 
45 #include "G4INCLSurfaceAvatar.hh"
46 #include "G4INCLRandom.hh"
49 #include "G4INCLClustering.hh"
50 #include <sstream>
51 #include <string>
52 
53 namespace G4INCL {
54 
56  :IAvatar(time), theParticle(aParticle), theNucleus(n),
57  particlePIn(0.),
58  particlePOut(0.),
59  particleTOut(0.),
60  TMinusV(0.),
61  TMinusV2(0.),
62  particleMass(0.),
63  sinIncidentAngle(0.),
64  cosIncidentAngle(0.),
65  sinRefractionAngle(0.),
66  cosRefractionAngle(0.),
67  refractionIndexRatio(0.),
68  internalReflection(false)
69  {
71  }
72 
74 
75  }
76 
78  if(theParticle->isTargetSpectator()) {
79  INCL_DEBUG("Particle " << theParticle->getID() << " is a spectator, reflection" << '\n');
80  return new ReflectionChannel(theNucleus, theParticle);
81  }
82 
83  // We forbid transmission of resonances below the Fermi energy. Emitting a
84  // delta particle below Tf can lead to negative excitation energies, since
85  // CDPP assumes that particles stay in the Fermi sea.
86  if(theParticle->isResonance()) {
87  const G4double theFermiEnergy = theNucleus->getPotential()->getFermiEnergy(theParticle);
88  if(theParticle->getKineticEnergy()<theFermiEnergy) {
89  INCL_DEBUG("Particle " << theParticle->getID() << " is a resonance below Tf, reflection" << '\n'
90  << " Tf=" << theFermiEnergy << ", EKin=" << theParticle->getKineticEnergy() << '\n');
91  return new ReflectionChannel(theNucleus, theParticle);
92  }
93  }
94 
95  // Don't try to make a cluster if the leading particle is too slow
96  const G4double transmissionProbability = getTransmissionProbability(theParticle);
97  const G4double TOut = TMinusV;
98  const G4double kOut = particlePOut;
99  const G4double cosR = cosRefractionAngle;
100 
101  INCL_DEBUG("Transmission probability for particle " << theParticle->getID() << " = " << transmissionProbability << '\n');
102  /* Don't attempt to construct clusters when a projectile spectator is
103  * trying to escape during a nucleus-nucleus collision. The idea behind
104  * this is that projectile spectators will later be collected in the
105  * projectile remnant, and trying to clusterise them somewhat feels like
106  * G4double counting. Moreover, applying the clustering algorithm on escaping
107  * projectile spectators makes the code *really* slow if the projectile is
108  * large.
109  */
110  if(theParticle->isNucleon()
111  && (!theParticle->isProjectileSpectator() || !theNucleus->isNucleusNucleusCollision())
112  && transmissionProbability>1.E-4) {
113  Cluster *candidateCluster = 0;
114 
115  candidateCluster = Clustering::getCluster(theNucleus, theParticle);
116  if(candidateCluster != 0 &&
117  Clustering::clusterCanEscape(theNucleus, candidateCluster)) {
118 
119  INCL_DEBUG("Cluster algorithm succeded. Candidate cluster:" << '\n' << candidateCluster->print() << '\n');
120 
121  // Check if the cluster can penetrate the Coulomb barrier
122  const G4double clusterTransmissionProbability = getTransmissionProbability(candidateCluster);
123  const G4double x = Random::shoot();
124 
125  INCL_DEBUG("Transmission probability for cluster " << candidateCluster->getID() << " = " << clusterTransmissionProbability << '\n');
126 
127  if (x <= clusterTransmissionProbability) {
128  theNucleus->getStore()->getBook().incrementEmittedClusters();
129  INCL_DEBUG("Cluster " << candidateCluster->getID() << " passes the Coulomb barrier, transmitting." << '\n');
130  return new TransmissionChannel(theNucleus, candidateCluster);
131  } else {
132  INCL_DEBUG("Cluster " << candidateCluster->getID() << " does not pass the Coulomb barrier. Falling back to transmission of the leading particle." << '\n');
133  delete candidateCluster;
134  }
135  } else {
136  delete candidateCluster;
137  }
138  }
139 
140  // If we haven't transmitted a cluster (maybe cluster feature was
141  // disabled or maybe we just can't produce an acceptable cluster):
142 
143  // Always transmit projectile spectators if no cluster was formed and if
144  // transmission is energetically allowed
145  if(theParticle->isProjectileSpectator() && transmissionProbability>0.) {
146  INCL_DEBUG("Particle " << theParticle->getID() << " is a projectile spectator, transmission" << '\n');
147  return new TransmissionChannel(theNucleus, theParticle, TOut);
148  }
149 
150  // Transmit or reflect depending on the transmission probability
151  const G4double x = Random::shoot();
152 
153  if(x <= transmissionProbability) { // Transmission
154  INCL_DEBUG("Particle " << theParticle->getID() << " passes the Coulomb barrier, transmitting." << '\n');
155  if(theNucleus->getStore()->getConfig()->getRefraction()) {
156  return new TransmissionChannel(theNucleus, theParticle, kOut, cosR);
157  } else {
158  return new TransmissionChannel(theNucleus, theParticle, TOut);
159  }
160  } else { // Reflection
161  INCL_DEBUG("Particle " << theParticle->getID() << " does not pass the Coulomb barrier, reflection." << '\n');
162  return new ReflectionChannel(theNucleus, theParticle);
163  }
164  }
165 
167  getChannel()->fillFinalState(fs);
168  }
169 
171 
173  ParticleList const &outgoing = fs->getOutgoingParticles();
174  if(!outgoing.empty()) { // Transmission
175 // assert(outgoing.size()==1);
176  Particle *out = outgoing.front();
177  out->rpCorrelate();
178  if(out->isCluster()) {
179  Cluster *clusterOut = dynamic_cast<Cluster*>(out);
180  ParticleList const &components = clusterOut->getParticles();
181  for(ParticleIter i=components.begin(), e=components.end(); i!=e; ++i) {
182  if(!(*i)->isTargetSpectator())
183  theNucleus->getStore()->getBook().decrementCascading();
184  }
185  } else if(!theParticle->isTargetSpectator()) {
186 // assert(out==theParticle);
187  theNucleus->getStore()->getBook().decrementCascading();
188  }
189  }
190  }
191 
192  std::string SurfaceAvatar::dump() const {
193  std::stringstream ss;
194  ss << "(avatar " << theTime << " 'reflection" << '\n'
195  << "(list " << '\n'
196  << theParticle->dump()
197  << "))" << '\n';
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  *(Math::arcCos(px)-px*std::sqrt(1.-px*px));
264  INCL_DEBUG("Coulomb barrier, logCoulombTransmission=" << logCoulombTransmission << '\n');
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 << '\n');
294  }
295 }
G4int getA() const
Returns the baryon number.
G4double arcCos(const G4double x)
Calculates arccos with some tolerance on illegal arguments.
G4double getCosRPAngle() const
Get the cosine of the angle between position and momentum.
G4bool isResonance() const
Is it a resonance?
G4double getMass() const
Get the cached particle mass.
ParticleList const & getParticles() const
G4bool isNucleusNucleusCollision() const
Is it a nucleus-nucleus collision?
G4bool isTargetSpectator() const
const G4INCL::ThreeVector & getMomentum() const
Config const * getConfig()
Definition: G4INCLStore.hh:273
Store * getStore() const
tuple x
Definition: test.py:50
virtual void postInteraction(FinalState *)
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.
G4double getFermiEnergy(const Particle *const p) const
Return the Fermi energy for a particle.
virtual void fillFinalState(FinalState *fs)=0
G4bool clusterCanEscape(Nucleus const *const n, Cluster const *const c)
Determine whether the cluster can escape or not.
void rpCorrelate()
Make the particle follow a strict r-p correlation.
ParticleList const & getOutgoingParticles() const
void setType(AvatarType t)
SurfaceAvatar(G4INCL::Particle *aParticle, G4double time, G4INCL::Nucleus *aNucleus)
Book & getBook()
Definition: G4INCLStore.hh:259
std::string print() const
G4int getZ() const
Returns the charge number.
void incrementEmittedClusters()
Definition: G4INCLBook.hh:79
const G4int n
Cluster * getCluster(Nucleus *n, Particle *p)
Call the clustering algorithm.
NuclearPotential::INuclearPotential const * getPotential() const
Getter for thePotential.
void decrementCascading()
Definition: G4INCLBook.hh:78
Static class for cluster formation.
G4double getTransmissionBarrier(Particle const *const p)
Get the transmission barrier.
G4bool isProjectileSpectator() const
std::string dump() const
G4bool isNucleon() const
G4bool isCluster() const
G4double shoot()
Definition: G4INCLRandom.cc:93
G4double getKineticEnergy() const
Get the particle kinetic energy.
double G4double
Definition: G4Types.hh:76
void fillFinalState(FinalState *fs)
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.