Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/hadronic/models/inclxx/incl_physics/src/G4INCLCoulombNonRelativistic.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

  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 /** \file G4INCLCoulombNonRelativistic.cc
 39  * \brief Class for non-relativistic Coulomb distortion.
 40  *
 41  * \date 14 February 2011
 42  * \author Davide Mancusi
 43  */
 44 
 45 #include "G4INCLCoulombNonRelativistic.hh"
 46 #include "G4INCLGlobals.hh"
 47 
 48 namespace G4INCL {
 49 
 50   ParticleEntryAvatar *CoulombNonRelativistic::bringToSurface(Particle * const p, Nucleus * const n) const {
 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
 60     return theCoulombNoneSlave.bringToSurface(p,n);
 61   }
 62 
 63   IAvatarList CoulombNonRelativistic::bringToSurface(Cluster * const c, Nucleus * const n) const {
 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
 75     return theCoulombNoneSlave.bringToSurface(c,n);
 76   }
 77 
 78   void CoulombNonRelativistic::distortOut(ParticleList const &pL,
 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 
124   G4double CoulombNonRelativistic::maxImpactParameter(ParticleSpecies const &p, const G4double kinE,
125                                                     Nucleus const * const n) const {
126     const G4double theMinimumDistance = minimumDistance(p, kinE, n);
127     G4double rMax = n->getUniverseRadius();
128     if(p.theType == Composite)
129       rMax +=  2.*ParticleTable::getLargestNuclearRadius(p.theA, p.theZ);
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 
137   G4bool CoulombNonRelativistic::coulombDeviation(Particle * const p, Nucleus const * const n) const {
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 
153     const G4double radius = getCoulombRadius(p->getSpecies(), n);
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 
196   G4double CoulombNonRelativistic::getCoulombRadius(ParticleSpecies const &p, Nucleus const * const n) const {
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.) {
223         radius = ParticleTable::getLargestNuclearRadius(Ap,Zp) + ParticleTable::getLargestNuclearRadius(At, Zt);
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 }
235