Geant4  v4-10.4-release
 모두 클래스 네임스페이스들 파일들 함수 변수 타입정의 열거형 타입 열거형 멤버 Friends 매크로 그룹들 페이지들
G4INCLCoulombNonRelativistic.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 
46 #include "G4INCLGlobals.hh"
47 
48 namespace G4INCL {
49 
51  // No distortion for neutral particles
52  if(p->getZ()!=0) {
53  const G4bool success = coulombDeviation(p, n);
54  if(!success) // transparent
55  return NULL;
56  }
57 
58  // Rely on the CoulombNone slave to compute the straight-line intersection
59  // and actually bring the particle to the surface of the nucleus
61  }
62 
64  // Neutral clusters?!
65 // assert(c->getZ()>0);
66 
67  // Perform the actual Coulomb deviation
68  const G4bool success = coulombDeviation(c, n);
69  if(!success) {
70  return IAvatarList();
71  }
72 
73  // Rely on the CoulombNone slave to compute the straight-line intersection
74  // and actually bring the particle to the surface of the nucleus
76  }
77 
79  Nucleus const * const nucleus) const {
80 
81  for(ParticleIter particle=pL.begin(), e=pL.end(); particle!=e; ++particle) {
82 
83  const G4int Z = (*particle)->getZ();
84  if(Z == 0) continue;
85 
86  const G4double tcos=1.-0.000001;
87 
88  const G4double et1 = PhysicalConstants::eSquared * nucleus->getZ();
89  const G4double transmissionRadius =
90  nucleus->getDensity()->getTransmissionRadius(*particle);
91 
92  const ThreeVector position = (*particle)->getPosition();
93  ThreeVector momentum = (*particle)->getMomentum();
94  const G4double r = position.mag();
95  const G4double p = momentum.mag();
96  const G4double cosTheta = position.dot(momentum)/(r*p);
97  if(cosTheta < 0.999999) {
98  const G4double sinTheta = std::sqrt(1.-cosTheta*cosTheta);
99  const G4double eta = et1 * Z / (*particle)->getKineticEnergy();
100  if(eta > transmissionRadius-0.0001) {
101  // If below the Coulomb barrier, radial emission:
102  momentum = position * (p/r);
103  (*particle)->setMomentum(momentum);
104  } else {
105  const G4double b0 = 0.5 * (eta + std::sqrt(eta*eta +
106  4. * std::pow(transmissionRadius*sinTheta,2)
107  * (1.-eta/transmissionRadius)));
108  const G4double bInf = std::sqrt(b0*(b0-eta));
109  const G4double thr = std::atan(eta/(2.*bInf));
110  G4double uTemp = (1.-b0/transmissionRadius) * std::sin(thr) +
111  b0/transmissionRadius;
112  if(uTemp>tcos) uTemp=tcos;
113  const G4double thd = Math::arcCos(cosTheta)-Math::piOverTwo + thr +
114  Math::arcCos(uTemp);
115  const G4double c1 = std::sin(thd)*cosTheta/sinTheta + std::cos(thd);
116  const G4double c2 = -p*std::sin(thd)/(r*sinTheta);
117  const ThreeVector newMomentum = momentum*c1 + position*c2;
118  (*particle)->setMomentum(newMomentum);
119  }
120  }
121  }
122  }
123 
125  Nucleus const * const n) const {
126  const G4double theMinimumDistance = minimumDistance(p, kinE, n);
127  G4double rMax = n->getUniverseRadius();
128  if(p.theType == Composite)
130  const G4double theMaxImpactParameterSquared = rMax*(rMax-theMinimumDistance);
131  if(theMaxImpactParameterSquared<=0.)
132  return 0.;
133  const G4double theMaxImpactParameter = std::sqrt(theMaxImpactParameterSquared);
134  return theMaxImpactParameter;
135  }
136 
138  // Determine the rotation angle and the new impact parameter
139  ThreeVector positionTransverse = p->getTransversePosition();
140  const G4double impactParameterSquared = positionTransverse.mag2();
141  const G4double impactParameter = std::sqrt(impactParameterSquared);
142 
143  // Some useful variables
144  const G4double theMinimumDistance = minimumDistance(p, n);
145  // deltaTheta2 = (pi - Rutherford scattering angle)/2
146  G4double deltaTheta2 = std::atan(2.*impactParameter/theMinimumDistance);
147  if(deltaTheta2<0.)
148  deltaTheta2 += Math::pi;
149  const G4double eccentricity = 1./std::cos(deltaTheta2);
150 
151  G4double newImpactParameter, alpha; // Parameters that must be determined by the deviation
152 
154  const G4double impactParameterTangentSquared = radius*(radius-theMinimumDistance);
155  if(impactParameterSquared >= impactParameterTangentSquared) {
156  // The particle trajectory misses the Coulomb sphere
157  // In this case the new impact parameter is the minimum distance of
158  // approach of the hyperbola
159 // assert(std::abs(1. + 2.*impactParameter*impactParameter/(radius*theMinimumDistance))>=eccentricity);
160  newImpactParameter = 0.5 * theMinimumDistance * (1.+eccentricity); // the minimum distance of approach
161  alpha = Math::piOverTwo - deltaTheta2; // half the Rutherford scattering angle
162  } else {
163  // The particle trajectory intersects the Coulomb sphere
164 
165  // Compute the entrance angle
166  const G4double argument = -(1. + 2.*impactParameter*impactParameter/(radius*theMinimumDistance))
167  / eccentricity;
168  const G4double thetaIn = Math::twoPi - Math::arcCos(argument) - deltaTheta2;
169 
170  // Velocity angle at the entrance point
171  alpha = std::atan((1+std::cos(thetaIn))
172  / (std::sqrt(eccentricity*eccentricity-1.) - std::sin(thetaIn)))
173  * Math::sign(theMinimumDistance);
174  // New impact parameter
175  newImpactParameter = radius * std::sin(thetaIn - alpha);
176  }
177 
178  // Modify the impact parameter of the particle
179  positionTransverse *= newImpactParameter/positionTransverse.mag();
180  const ThreeVector theNewPosition = p->getLongitudinalPosition() + positionTransverse;
181  p->setPosition(theNewPosition);
182 
183  // Determine the rotation axis for the incoming particle
184  const ThreeVector &momentum = p->getMomentum();
185  ThreeVector rotationAxis = momentum.vector(positionTransverse);
186  const G4double axisLength = rotationAxis.mag();
187  // Apply the rotation
188  if(axisLength>1E-20) {
189  rotationAxis /= axisLength;
190  p->rotatePositionAndMomentum(alpha, rotationAxis);
191  }
192 
193  return true;
194  }
195 
197  if(p.theType == Composite) {
198  const G4int Zp = p.theZ;
199  const G4int Ap = p.theA;
200  const G4int Zt = n->getZ();
201  const G4int At = n->getA();
202  G4double barr, radius = 0.;
203  if(Zp==1 && Ap==2) { // d
204  barr = 0.2565*Math::pow23((G4double)At)-0.78;
205  radius = PhysicalConstants::eSquared*Zp*Zt/barr - 2.5;
206  } else if(Zp==1 && Ap==3) { // t
207  barr = 0.5*(0.5009*Math::pow23((G4double)At)-1.16);
208  radius = PhysicalConstants::eSquared*Zt/barr - 0.5;
209  } else if(Zp==2) { // alpha, He3
210  barr = 0.5939*Math::pow23((G4double)At)-1.64;
211  radius = PhysicalConstants::eSquared*Zp*Zt/barr - 0.5;
212  } else if(Zp>2) {
213  // Coulomb radius from the Shen model
214  const G4double Ap13 = Math::pow13((G4double)Ap);
215  const G4double At13 = Math::pow13((G4double)At);
216  const G4double rp = 1.12*Ap13 - 0.94/Ap13;
217  const G4double rt = 1.12*At13 - 0.94/At13;
218  const G4double someRadius = rp+rt+3.2;
219  const G4double theShenBarrier = PhysicalConstants::eSquared*Zp*Zt/someRadius - rt*rp/(rt+rp);
220  radius = PhysicalConstants::eSquared*Zp*Zt/theShenBarrier;
221  }
222  if(radius<=0.) {
224  INCL_ERROR("Negative Coulomb radius! Using the sum of nuclear radii = " << radius << '\n');
225  }
226  INCL_DEBUG("Coulomb radius for particle "
227  << ParticleTable::getShortName(p) << " in nucleus A=" << At <<
228  ", Z=" << Zt << ": " << radius << '\n');
229  return radius;
230  } else
231  return n->getUniverseRadius();
232  }
233 
234 }
G4double getLargestNuclearRadius(const G4int A, const G4int Z)
G4double getCoulombRadius(ParticleSpecies const &p, Nucleus const *const n) const
Get the Coulomb radius for a given particle.
ParticleList::const_iterator ParticleIter
const G4double piOverTwo
virtual G4INCL::ParticleSpecies getSpecies() const
Get the particle species.
G4double getTransmissionRadius(Particle const *const p) const
The radius used for calculating the transmission coefficient.
ThreeVector vector(const ThreeVector &v) const
const G4INCL::ThreeVector & getMomentum() const
Class for non-relativistic Coulomb distortion.
G4double dot(const ThreeVector &v) const
const char * p
Definition: xmltok.h:285
void distortOut(ParticleList const &pL, Nucleus const *const n) const
Modify the momenta of the outgoing particles.
NuclearDensity const * getDensity() const
Getter for theDensity.
#define INCL_DEBUG(x)
Float_t Z
double G4double
Definition: G4Types.hh:76
bool G4bool
Definition: G4Types.hh:79
G4double minimumDistance(ParticleSpecies const &p, const G4double kineticEnergy, Nucleus const *const n) const
Return the minimum distance of approach in a head-on collision (b=0).
TCanvas * c2
Definition: plot_hist.C:75
virtual void rotatePositionAndMomentum(const G4double angle, const ThreeVector &axis)
Rotate the particle position and momentum.
virtual void setPosition(const G4INCL::ThreeVector &position)
#define INCL_ERROR(x)
G4double pow13(G4double x)
static const G4double alpha
G4bool coulombDeviation(Particle *const p, Nucleus const *const n) const
Perform Coulomb deviation.
G4double arcCos(const G4double x)
Calculates arccos with some tolerance on illegal arguments.
G4int getZ() const
Returns the charge number.
G4double getUniverseRadius() const
Getter for theUniverseRadius.
G4int getA() const
Returns the baryon number.
Double_t radius
const G4double pi
G4double maxImpactParameter(ParticleSpecies const &p, const G4double kinE, Nucleus const *const n) const
Return the maximum impact parameter for Coulomb-distorted trajectories.
CoulombNone theCoulombNoneSlave
Internal CoulombNone slave to generate the avatars.
G4int sign(const T t)
const G4double twoPi
ParticleEntryAvatar * bringToSurface(Particle *const p, Nucleus *const n) const
Modify the momentum of the particle and position it on the surface of the nucleus.
int G4int
Definition: G4Types.hh:78
ThreeVector getLongitudinalPosition() const
Longitudinal component of the position w.r.t. the momentum.
std::string getShortName(const ParticleType t)
Get the short INCL name of the particle.
G4double pow23(G4double x)
G4double mag2() const
G4double mag() const
const G4double eSquared
Coulomb conversion factor [MeV*fm].
ParticleEntryAvatar * bringToSurface(Particle *const p, Nucleus *const n) const
Position the particle on the surface of the nucleus.
Char_t n[5]
ThreeVector getTransversePosition() const
Transverse component of the position w.r.t. the momentum.
UnorderedVector< IAvatar * > IAvatarList