00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033 #define INCLXX_IN_GEANT4_MODE 1
00034
00035 #include "globals.hh"
00036
00048 #ifndef G4INCLINUCLEARPOTENTIAL_HH
00049 #define G4INCLINUCLEARPOTENTIAL_HH 1
00050
00051 #include "G4INCLParticle.hh"
00052 #include "G4INCLRandom.hh"
00053 #include "G4INCLDeuteronDensity.hh"
00054 #include <map>
00055
00056
00057 namespace G4INCL {
00058
00059 namespace NuclearPotential {
00060
00061 class INuclearPotential {
00062 public:
00063 INuclearPotential(const G4int A, const G4int Z, const G4bool pionPot) :
00064 theA(A),
00065 theZ(Z),
00066 pionPotential(pionPot)
00067 {
00068 if(pionPotential) {
00069 const G4double ZOverA = ((G4double) theZ) / ((G4double) theA);
00070
00071 const G4double r = 1.12*Math::pow13((G4double)theA);
00072
00073 const G4double xsi = 1. - 2.*ZOverA;
00074 const G4double vc = 1.25*PhysicalConstants::eSquared*theZ/r;
00075 vPiPlus = vPionDefault + 71.*xsi - vc;
00076 vPiZero = vPionDefault;
00077 vPiMinus = vPionDefault - 71.*xsi + vc;
00078 } else {
00079 vPiPlus = 0.0;
00080 vPiZero = 0.0;
00081 vPiMinus = 0.0;
00082 }
00083 }
00084
00085 virtual ~INuclearPotential() {}
00086
00088 G4bool hasPionPotential() { return pionPotential; }
00089
00090 virtual G4double computePotentialEnergy(const Particle * const p) const = 0;
00091
00097 inline G4double getFermiEnergy(const Particle * const p) const {
00098 std::map<ParticleType, G4double>::const_iterator i = fermiEnergy.find(p->getType());
00099
00100 return i->second;
00101 }
00102
00108 inline G4double getFermiEnergy(const ParticleType t) const {
00109 std::map<ParticleType, G4double>::const_iterator i = fermiEnergy.find(t);
00110
00111 return i->second;
00112 }
00113
00119 inline G4double getSeparationEnergy(const Particle * const p) const {
00120 std::map<ParticleType, G4double>::const_iterator i = separationEnergy.find(p->getType());
00121
00122 return i->second;
00123 }
00124
00130 inline G4double getSeparationEnergy(const ParticleType t) const {
00131 std::map<ParticleType, G4double>::const_iterator i = separationEnergy.find(t);
00132
00133 return i->second;
00134 }
00135
00141 inline G4double getFermiMomentum(const Particle * const p) const {
00142 if(p->isDelta()) {
00143 const G4double Tf = getFermiEnergy(p), m = p->getMass();
00144 return std::sqrt(Tf*(Tf+2.*m));
00145 } else {
00146 std::map<ParticleType, G4double>::const_iterator i = fermiMomentum.find(p->getType());
00147
00148 return i->second;
00149 }
00150 }
00151
00157 inline G4double getFermiMomentum(const ParticleType t) const {
00158
00159 std::map<ParticleType, G4double>::const_iterator i = fermiMomentum.find(t);
00160 return i->second;
00161 }
00162
00163 protected:
00165 G4double computePionPotentialEnergy(const Particle * const p) const {
00166
00167 if(pionPotential && !p->isOutOfWell()) {
00168 switch( p->getType() ) {
00169 case PiPlus:
00170 return vPiPlus;
00171 break;
00172 case PiZero:
00173 return vPiZero;
00174 break;
00175 case PiMinus:
00176 return vPiMinus;
00177 break;
00178 default:
00179 return 0.0;
00180 break;
00181 }
00182 }
00183 else
00184 return 0.0;
00185 }
00186
00187 protected:
00189 const G4int theA;
00191 const G4int theZ;
00192 private:
00193 const G4bool pionPotential;
00194 G4double vPiPlus, vPiZero, vPiMinus;
00195 static const G4double vPionDefault;
00196 protected:
00197
00198 std::map<ParticleType,G4double> fermiEnergy;
00199
00200 std::map<ParticleType,G4double> fermiMomentum;
00201
00202 std::map<ParticleType,G4double> separationEnergy;
00203
00204 };
00205
00206 }
00207
00208 }
00209
00210 #endif