| Geant4 Cross Reference |
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 G4INCLINuclearPotential.hh
39 * \brief Abstract interface to the nuclear potential.
40 *
41 * NuclearPotential-like classes should provide access to the value of the
42 * potential of a particle in a particular context. For example, an instance of
43 * a NuclearPotential class should be associated to every nucleus.
44 *
45 * \date 17 January 2011
46 * \author Davide Mancusi
47 */
48
49 #ifndef G4INCLINUCLEARPOTENTIAL_HH
50 #define G4INCLINUCLEARPOTENTIAL_HH 1
51
52 #include "G4INCLParticle.hh"
53 #include "G4INCLRandom.hh"
54 #include "G4INCLDeuteronDensity.hh"
55 #include <map>
56 // #include <cassert>
57
58 namespace G4INCL {
59
60 namespace NuclearPotential {
61 class INuclearPotential {
62 public:
63 INuclearPotential(const G4int A, const G4int Z, const G4bool pionPot) :
64 theA(A),
65 theZ(Z),
66 pionPotential(pionPot)
67 {
68 if(pionPotential) {
69 const G4double ZOverA = ((G4double) theZ) / ((G4double) theA);
70 // As in INCL4.6, use the r0*A^(1/3) formula to estimate vc
71 const G4double r = 1.12*Math::pow13((G4double)theA);
72
73 const G4double xsi = 1. - 2.*ZOverA;
74 const G4double vc = 1.25*PhysicalConstants::eSquared*theZ/r;
75 vPiPlus = vPionDefault + 71.*xsi - vc;
76 vPiZero = vPionDefault;
77 vPiMinus = vPionDefault - 71.*xsi + vc;
78 vKPlus = vKPlusDefault;
79 vKZero = vKPlusDefault + 10.; // Hypothesis to be check
80 vKMinus = vKMinusDefault;
81 vKZeroBar = vKMinusDefault - 10.; // Hypothesis to be check
82 } else {
83 vPiPlus = 0.0;
84 vPiZero = 0.0;
85 vPiMinus = 0.0;
86 vKPlus = 0.0;
87 vKZero = 0.0;
88 vKMinus = 0.0;
89 vKZeroBar = 0.0;
90 }
91 }
92
93 virtual ~INuclearPotential() {}
94
95 /// \brief Do we have a pion potential?
96 G4bool hasPionPotential() const { return pionPotential; }
97
98 virtual G4double computePotentialEnergy(const Particle * const p) const = 0;
99
100 /** \brief Return the Fermi energy for a particle.
101 *
102 * \param p pointer to a Particle
103 * \return Fermi energy for that particle type
104 **/
105 inline G4double getFermiEnergy(const Particle * const p) const {
106 std::map<ParticleType, G4double>::const_iterator i = fermiEnergy.find(p->getType());
107 // assert(i!=fermiEnergy.end());
108 return i->second;
109 }
110
111 /** \brief Return the Fermi energy for a particle type.
112 *
113 * \param t particle type
114 * \return Fermi energy for that particle type
115 **/
116 inline G4double getFermiEnergy(const ParticleType t) const {
117 std::map<ParticleType, G4double>::const_iterator i = fermiEnergy.find(t);
118 // assert(i!=fermiEnergy.end());
119 return i->second;
120 }
121
122 /** \brief Return the separation energy for a particle.
123 *
124 * \param p pointer to a Particle
125 * \return separation energy for that particle type
126 **/
127 inline G4double getSeparationEnergy(const Particle * const p) const {
128 std::map<ParticleType, G4double>::const_iterator i = separationEnergy.find(p->getType());
129 // assert(i!=separationEnergy.end());
130 return i->second;
131 }
132
133 /** \brief Return the separation energy for a particle type.
134 *
135 * \param t particle type
136 * \return separation energy for that particle type
137 **/
138 inline G4double getSeparationEnergy(const ParticleType t) const {
139 std::map<ParticleType, G4double>::const_iterator i = separationEnergy.find(t);
140 // assert(i!=separationEnergy.end());
141 return i->second;
142 }
143
144 /** \brief Return the Fermi momentum for a particle.
145 *
146 * \param p pointer to a Particle
147 * \return Fermi momentum for that particle type
148 **/
149 inline G4double getFermiMomentum(const Particle * const p) const {
150 if(p->isDelta()) {
151 const G4double Tf = getFermiEnergy(p), mass = p->getMass();
152 return std::sqrt(Tf*(Tf+2.*mass));
153 } else {
154 std::map<ParticleType, G4double>::const_iterator i = fermiMomentum.find(p->getType());
155 // assert(i!=fermiMomentum.end());
156 return i->second;
157 }
158 }
159
160 /** \brief Return the Fermi momentum for a particle type.
161 *
162 * \param t particle type
163 * \return Fermi momentum for that particle type
164 **/
165 inline G4double getFermiMomentum(const ParticleType t) const {
166 // assert(t!=DeltaPlusPlus && t!=DeltaPlus && t!=DeltaZero && t!=DeltaMinus);
167 std::map<ParticleType, G4double>::const_iterator i = fermiMomentum.find(t);
168 return i->second;
169 }
170
171 protected:
172 /// \brief Compute the potential energy for the given pion.
173 G4double computePionPotentialEnergy(const Particle * const p) const {
174 // assert(p->getType()==PiPlus || p->getType()==PiZero || p->getType()==PiMinus);
175 if(pionPotential && !p->isOutOfWell()) {
176 switch( p->getType() ) {
177 case PiPlus:
178 return vPiPlus;
179 break;
180 case PiZero:
181 return vPiZero;
182 break;
183 case PiMinus:
184 return vPiMinus;
185 break;
186 default: // Pion potential is defined and non-zero only for pions
187 return 0.0;
188 break;
189 }
190 }
191 else
192 return 0.0;
193 }
194
195 protected:
196 /// \brief Compute the potential energy for the given kaon.
197 G4double computeKaonPotentialEnergy(const Particle * const p) const {
198 // assert(p->getType()==KPlus || p->getType()==KZero || p->getType()==KZeroBar || p->getType()==KMinus|| p->getType()==KShort|| p->getType()==KLong);
199 if(pionPotential && !p->isOutOfWell()) { // if pionPotental false -> kaonPotential false
200 switch( p->getType() ) {
201 case KPlus:
202 return vKPlus;
203 break;
204 case KZero:
205 return vKZero;
206 break;
207 case KZeroBar:
208 return vKZeroBar;
209 break;
210 case KShort:
211 case KLong:
212 return 0.0; // Should never be in the nucleus
213 break;
214 case KMinus:
215 return vKMinus;
216 break;
217 default:
218 return 0.0;
219 break;
220 }
221 }
222 else
223 return 0.0;
224 }
225
226 protected:
227 /// \brief Compute the potential energy for the given pion resonances (Eta, Omega and EtaPrime and Gamma also).
228 G4double computePionResonancePotentialEnergy(const Particle * const p) const {
229 // assert(p->getType()==Eta || p->getType()==Omega || p->getType()==EtaPrime || p->getType()==Photon);
230 if(pionPotential && !p->isOutOfWell()) {
231 switch( p->getType() ) {
232 case Eta:
233 //jcd return vPiZero;
234 //jcd return vPiZero*1.5;
235 return 0.0; // (JCD: seems to give better results)
236 break;
237 case Omega:
238 return 15.0; // S.Friedrich et al., Physics Letters B736(2014)26-32. (V. Metag in Hyperfine Interact (2015) 234:25-31 gives 29 MeV)
239 break;
240 case EtaPrime:
241 return 37.0; // V. Metag in Hyperfine Interact (2015) 234:25-31
242 break;
243 case Photon:
244 return 0.0;
245 break;
246 default:
247 return 0.0;
248 break;
249 }
250 }
251 else
252 return 0.0;
253 }
254
255 protected:
256 /// \brief The mass number of the nucleus
257 const G4int theA;
258 /// \brief The charge number of the nucleus
259 const G4int theZ;
260 private:
261 const G4bool pionPotential;
262 G4double vPiPlus, vPiZero, vPiMinus;
263 static const G4double vPionDefault;
264 G4double vKPlus, vKZero, vKZeroBar, vKMinus;
265 static const G4double vKPlusDefault;
266 static const G4double vKMinusDefault;
267 protected:
268 /* \brief map of Fermi energies per particle type */
269 std::map<ParticleType,G4double> fermiEnergy;
270 /* \brief map of Fermi momenta per particle type */
271 std::map<ParticleType,G4double> fermiMomentum;
272 /* \brief map of separation energies per particle type */
273 std::map<ParticleType,G4double> separationEnergy;
274
275 };
276
277
278
279 /** \brief Create an INuclearPotential object
280 *
281 * This is the method that should be used to instantiate objects derived
282 * from INuclearPotential. It uses a caching mechanism to minimise
283 * thrashing and speed up the code.
284 *
285 * \param type the type of the potential to be created
286 * \param theA mass number of the nucleus
287 * \param theZ charge number of the nucleus
288 * \param pionPotential whether pions should also feel the potential
289 * \return a pointer to the nuclear potential
290 */
291 INuclearPotential const *createPotential(const PotentialType type, const G4int theA, const G4int theZ, const G4bool pionPotential);
292
293 /// \brief Clear the INuclearPotential cache
294 void clearCache();
295
296 }
297
298 }
299
300 #endif /* G4INCLINUCLEARPOTENTIAL_HH_ */
301