ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4INCLINuclearPotential.hh
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4INCLINuclearPotential.hh
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 
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 {
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;
75  vPiPlus = vPionDefault + 71.*xsi - vc;
77  vPiMinus = vPionDefault - 71.*xsi + vc;
79  vKZero = vKPlusDefault + 10.; // Hypothesis to be check
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 
97 
98  virtual G4double computePotentialEnergy(const Particle * const p) const = 0;
99 
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 
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 
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 
139  std::map<ParticleType, G4double>::const_iterator i = separationEnergy.find(t);
140 // assert(i!=separationEnergy.end());
141  return i->second;
142  }
143 
149  inline G4double getFermiMomentum(const Particle * const p) const {
150  if(p->isDelta()) {
151  const G4double Tf = getFermiEnergy(p), m = p->getMass();
152  return std::sqrt(Tf*(Tf+2.*m));
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 
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:
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:
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:
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:
257  const G4int theA;
259  const G4int theZ;
260  private:
263  static const G4double vPionDefault;
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 
291  INuclearPotential const *createPotential(const PotentialType type, const G4int theA, const G4int theZ, const G4bool pionPotential);
292 
294  void clearCache();
295 
296  }
297 
298 }
299 
300 #endif /* G4INCLINUCLEARPOTENTIAL_HH_ */