ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4INCLNuclearPotentialConstant.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4INCLNuclearPotentialConstant.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 
47 #include "G4INCLNuclearPotentialConstant.hh"
48 #include "G4INCLParticleTable.hh"
49 
50 namespace G4INCL {
51 
52  namespace NuclearPotential {
53 
54  // Constructors
55  NuclearPotentialConstant::NuclearPotentialConstant(const G4int A, const G4int Z, const G4bool aPionPotential)
56  : INuclearPotential(A, Z, aPionPotential)
57  {
58  initialize();
59  }
60 
61  // Destructor
62  NuclearPotentialConstant::~NuclearPotentialConstant() {
63  }
64 
65  void NuclearPotentialConstant::initialize() {
66  const G4double mp = ParticleTable::getINCLMass(Proton);
67  const G4double mn = ParticleTable::getINCLMass(Neutron);
68 
69  const G4double theFermiMomentum = ParticleTable::getFermiMomentum(theA,theZ);
70 
71  fermiMomentum[Proton] = theFermiMomentum;
72  const G4double theProtonFermiEnergy = std::sqrt(theFermiMomentum*theFermiMomentum + mp*mp) - mp;
73  fermiEnergy[Proton] = theProtonFermiEnergy;
74 
75  fermiMomentum[Neutron] = theFermiMomentum;
76  const G4double theNeutronFermiEnergy = std::sqrt(theFermiMomentum*theFermiMomentum + mn*mn) - mn;
77  fermiEnergy[Neutron] = theNeutronFermiEnergy;
78 
79  fermiEnergy[DeltaPlusPlus] = fermiEnergy.find(Proton)->second;
80  fermiEnergy[DeltaPlus] = fermiEnergy.find(Proton)->second;
81  fermiEnergy[DeltaZero] = fermiEnergy.find(Neutron)->second;
82  fermiEnergy[DeltaMinus] = fermiEnergy.find(Neutron)->second;
83 
84  fermiEnergy[SigmaPlus] = fermiEnergy.find(Proton)->second;
85  fermiEnergy[SigmaZero] = fermiEnergy.find(Proton)->second;
86  fermiEnergy[SigmaMinus] = fermiEnergy.find(Proton)->second;
87 
88  fermiEnergy[Lambda] = fermiEnergy.find(Neutron)->second;
89 
90 
91  const G4double theAverageSeparationEnergy = 0.5*(ParticleTable::getSeparationEnergy(Proton,theA,theZ)+ParticleTable::getSeparationEnergy(Neutron,theA,theZ));
92  separationEnergy[Proton] = theAverageSeparationEnergy;
93  separationEnergy[Neutron] = theAverageSeparationEnergy;
94 
95  // Use separation energies from the ParticleTable
96  vNucleon = 0.5*(theProtonFermiEnergy + theNeutronFermiEnergy) + theAverageSeparationEnergy;
97  vDelta = vNucleon;
98  vSigma = -16.; // Caution: repulsive potential for Sigmas
99  vLambda = 28.;
100  separationEnergy[DeltaPlusPlus] = vDelta - fermiEnergy.find(DeltaPlusPlus)->second;
101  separationEnergy[DeltaPlus] = vDelta - fermiEnergy.find(DeltaPlus)->second;
102  separationEnergy[DeltaZero] = vDelta - fermiEnergy.find(DeltaZero)->second;
103  separationEnergy[DeltaMinus] = vDelta - fermiEnergy.find(DeltaMinus)->second;
104 
105  separationEnergy[PiPlus] = 0.;
106  separationEnergy[PiZero] = 0.;
107  separationEnergy[PiMinus] = 0.;
108 
109  separationEnergy[Eta] = 0.;
110  separationEnergy[Omega] = 0.;
111  separationEnergy[EtaPrime] = 0.;
112  separationEnergy[Photon] = 0.;
113 
114  separationEnergy[Lambda] = 0.;
115  separationEnergy[SigmaPlus] = 0.;
116  separationEnergy[SigmaZero] = 0.;
117  separationEnergy[SigmaMinus] = 0.;
118  separationEnergy[KPlus] = 0.;
119  separationEnergy[KZero] = 0.;
120  separationEnergy[KZeroBar] = 0.;
121  separationEnergy[KMinus] = 0.;
122  separationEnergy[KShort] = 0.;
123  separationEnergy[KLong] = 0.;
124 
125  INCL_DEBUG("Table of separation energies [MeV] for A=" << theA << ", Z=" << theZ << ":" << '\n'
126  << " proton: " << separationEnergy[Proton] << '\n'
127  << " neutron: " << separationEnergy[Neutron] << '\n'
128  << " delta++: " << separationEnergy[DeltaPlusPlus] << '\n'
129  << " delta+: " << separationEnergy[DeltaPlus] << '\n'
130  << " delta0: " << separationEnergy[DeltaZero] << '\n'
131  << " delta-: " << separationEnergy[DeltaMinus] << '\n'
132  << " pi+: " << separationEnergy[PiPlus] << '\n'
133  << " pi0: " << separationEnergy[PiZero] << '\n'
134  << " pi-: " << separationEnergy[PiMinus] << '\n'
135  << " eta: " << separationEnergy[Eta] << '\n'
136  << " omega: " << separationEnergy[Omega] << '\n'
137  << " etaprime:" << separationEnergy[EtaPrime] << '\n'
138  << " photon: " << separationEnergy[Photon] << '\n'
139  << " lambda: " << separationEnergy[Lambda] << '\n'
140  << " sigmaplus: " << separationEnergy[SigmaPlus] << '\n'
141  << " sigmazero: " << separationEnergy[SigmaZero] << '\n'
142  << " sigmaminus: " << separationEnergy[SigmaMinus] << '\n'
143  << " kplus: " << separationEnergy[KPlus] << '\n'
144  << " kzero: " << separationEnergy[KZero] << '\n'
145  << " kzerobar: " << separationEnergy[KZeroBar] << '\n'
146  << " kminus: " << separationEnergy[KMinus] << '\n'
147  << " kshort: " << separationEnergy[KShort] << '\n'
148  << " klong: " << separationEnergy[KLong] << '\n'
149  );
150 
151  INCL_DEBUG("Table of Fermi energies [MeV] for A=" << theA << ", Z=" << theZ << ":" << '\n'
152  << " proton: " << fermiEnergy[Proton] << '\n'
153  << " neutron: " << fermiEnergy[Neutron] << '\n'
154  << " delta++: " << fermiEnergy[DeltaPlusPlus] << '\n'
155  << " delta+: " << fermiEnergy[DeltaPlus] << '\n'
156  << " delta0: " << fermiEnergy[DeltaZero] << '\n'
157  << " delta-: " << fermiEnergy[DeltaMinus] << '\n'
158  << " lambda: " << fermiEnergy[Lambda] << '\n'
159  << " sigmaplus: " << fermiEnergy[SigmaPlus] << '\n'
160  << " sigmazero: " << fermiEnergy[SigmaZero] << '\n'
161  << " sigmaminus: " << fermiEnergy[SigmaMinus] << '\n'
162  );
163 
164  INCL_DEBUG("Table of Fermi momenta [MeV/c] for A=" << theA << ", Z=" << theZ << ":" << '\n'
165  << " proton: " << fermiMomentum[Proton] << '\n'
166  << " neutron: " << fermiMomentum[Neutron] << '\n'
167  );
168  }
169 
170  G4double NuclearPotentialConstant::computePotentialEnergy(const Particle *particle) const {
171 
172  switch( particle->getType() )
173  {
174  case Proton:
175  case Neutron:
176  return vNucleon;
177  break;
178 
179  case PiPlus:
180  case PiZero:
181  case PiMinus:
182  return computePionPotentialEnergy(particle);
183  break;
184 
185  case Eta:
186  case Omega:
187  case EtaPrime:
188  return computePionResonancePotentialEnergy(particle);
189  break;
190 
191  case SigmaPlus:
192  case SigmaZero:
193  case SigmaMinus:
194  return vSigma;
195  break;
196  case Lambda:
197  return vLambda;
198  break;
199 
200  case KPlus:
201  case KZero:
202  case KZeroBar:
203  case KMinus:
204  return computeKaonPotentialEnergy(particle);
205  break;
206 
207  case Photon:
208  return 0.0;
209  break;
210 
211  case DeltaPlusPlus:
212  case DeltaPlus:
213  case DeltaZero:
214  case DeltaMinus:
215  return vDelta;
216  break;
217  case UnknownParticle:
218  INCL_ERROR("Trying to compute potential energy of an unknown particle.");
219  return 0.0;
220  break;
221  default:
222  INCL_ERROR("Trying to compute potential energy of a malformed particle.");
223  return 0.0;
224  break;
225  }
226  }
227 
228  }
229 }
230