ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4INCLCoulombNonRelativistic.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file 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 =
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 }