Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/hadronic/models/inclxx/incl_physics/src/G4INCLNuclearPotentialIsospin.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 G4INCLNuclearPotentialIsospin.cc
 39  * \brief Isospin-dependent nuclear potential.
 40  *
 41  * Provides an isospin-dependent nuclear potential.
 42  *
 43  * \date 28 February 2011
 44  * \author Davide Mancusi
 45  */
 46 
 47 #include "G4INCLNuclearPotentialIsospin.hh"
 48 #include "G4INCLNuclearPotentialConstant.hh"
 49 #include "G4INCLParticleTable.hh"
 50 #include "G4INCLGlobals.hh"
 51 
 52 namespace G4INCL {
 53 
 54   namespace NuclearPotential {
 55 
 56     // Constructors
 57     NuclearPotentialIsospin::NuclearPotentialIsospin(const G4int A, const G4int Z, const G4bool aPionPotential)
 58       : INuclearPotential(A, Z, aPionPotential)
 59     {
 60       initialize();
 61     }
 62 
 63     // Destructor
 64     NuclearPotentialIsospin::~NuclearPotentialIsospin() {}
 65 
 66     void NuclearPotentialIsospin::initialize() {
 67       const G4double ZOverA = ((G4double) theZ) / ((G4double) theA);
 68 
 69       const G4double mp = ParticleTable::getINCLMass(Proton);
 70       const G4double mn = ParticleTable::getINCLMass(Neutron);
 71       const G4double ml = ParticleTable::getINCLMass(Lambda);
 72 
 73       const G4double theFermiMomentum = ParticleTable::getFermiMomentum(theA,theZ);
 74 
 75       fermiMomentum[Proton] = theFermiMomentum * Math::pow13(2.*ZOverA);
 76       const G4double theProtonFermiEnergy = std::sqrt(fermiMomentum[Proton]*fermiMomentum[Proton] + mp*mp) - mp;
 77       fermiEnergy[Proton] = theProtonFermiEnergy;
 78       // Use separation energies from the ParticleTable
 79       const G4double theProtonSeparationEnergy = ParticleTable::getSeparationEnergy(Proton,theA,theZ);
 80       separationEnergy[Proton] = theProtonSeparationEnergy;
 81       vProton = theProtonFermiEnergy + theProtonSeparationEnergy;
 82 
 83       fermiMomentum[Neutron] = theFermiMomentum * Math::pow13(2.*(1.-ZOverA));
 84       const G4double theNeutronFermiEnergy = std::sqrt(fermiMomentum[Neutron]*fermiMomentum[Neutron] + mn*mn) - mn;
 85       fermiEnergy[Neutron] = theNeutronFermiEnergy;
 86       // Use separation energies from the ParticleTable
 87       const G4double theNeutronSeparationEnergy = ParticleTable::getSeparationEnergy(Neutron,theA,theZ);
 88       separationEnergy[Neutron] = theNeutronSeparationEnergy;
 89       vNeutron = theNeutronFermiEnergy + theNeutronSeparationEnergy;
 90 
 91       const G4double separationEnergyDeltaPlusPlus = 2.*theProtonSeparationEnergy - theNeutronSeparationEnergy;
 92       separationEnergy[DeltaPlusPlus] = separationEnergyDeltaPlusPlus;
 93       separationEnergy[DeltaPlus] = theProtonSeparationEnergy;
 94       separationEnergy[DeltaZero] = theNeutronSeparationEnergy;
 95       const G4double separationEnergyDeltaMinus = 2.*theNeutronSeparationEnergy - theProtonSeparationEnergy;
 96       separationEnergy[DeltaMinus] = separationEnergyDeltaMinus;
 97 
 98       const G4double tinyMargin = 1E-7;
 99       vDeltaPlus = vProton;
100       vDeltaZero = vNeutron;
101       vDeltaPlusPlus = std::max(separationEnergyDeltaPlusPlus + tinyMargin, 2.*vDeltaPlus - vDeltaZero);
102       vDeltaMinus = std::max(separationEnergyDeltaMinus + tinyMargin, 2.*vDeltaZero - vDeltaPlus);
103       
104       vSigmaMinus = -16.; // Repulsive potential, from Eur. Phys.J.A. (2016) 52:21
105       vSigmaZero = -16.;  // hypothesis: same potential for each sigma
106       vSigmaPlus = -16.;
107 
108       vLambda = 30.;
109       vantiProton = 100.;
110 
111       const G4double asy = (theA - 2.*theZ)/theA;
112       // Jose Luis Rodriguez-Sanchez et al., Rapid Communication PRC 98, 021602 (2018)
113       if      (asy > 0.236) vLambda = 40.91;
114       else if (asy > 0.133) vLambda = 56.549 - 678.73*asy + 4905.35*asy*asy - 9789.1*asy*asy*asy;
115          
116       const G4double theLambdaSeparationEnergy = ParticleTable::getSeparationEnergy(Lambda,theA,theZ);
117       const G4double theantiProtonSeparationEnergy = ParticleTable::getSeparationEnergy(antiProton,theA,theZ);
118 
119       separationEnergy[PiPlus] = theProtonSeparationEnergy - theNeutronSeparationEnergy;
120       separationEnergy[PiZero] = 0.;
121       separationEnergy[PiMinus] = theNeutronSeparationEnergy - theProtonSeparationEnergy;
122 
123       separationEnergy[Eta]      = 0.;
124       separationEnergy[Omega]    = 0.;
125       separationEnergy[EtaPrime] = 0.;
126       separationEnergy[Photon]   = 0.;
127       
128       separationEnergy[Lambda]    = theLambdaSeparationEnergy;
129       separationEnergy[SigmaPlus] = theProtonSeparationEnergy + theLambdaSeparationEnergy - theNeutronSeparationEnergy;
130       separationEnergy[SigmaZero] = theLambdaSeparationEnergy;
131       separationEnergy[SigmaMinus]  = theNeutronSeparationEnergy + theLambdaSeparationEnergy - theProtonSeparationEnergy;
132 
133       separationEnergy[KPlus]   = theProtonSeparationEnergy - theLambdaSeparationEnergy;
134       separationEnergy[KZero]   = (theNeutronSeparationEnergy - theLambdaSeparationEnergy);
135       separationEnergy[KZeroBar]  = (theLambdaSeparationEnergy - theNeutronSeparationEnergy);
136       separationEnergy[KMinus]    = 2.*theNeutronSeparationEnergy - theProtonSeparationEnergy-theLambdaSeparationEnergy;
137 
138       separationEnergy[KShort]    = (theNeutronSeparationEnergy - theLambdaSeparationEnergy);
139       separationEnergy[KLong]   = (theNeutronSeparationEnergy - theLambdaSeparationEnergy);
140 
141       separationEnergy[antiProton]    = theantiProtonSeparationEnergy;
142 
143       fermiEnergy[DeltaPlusPlus] = vDeltaPlusPlus - separationEnergy[DeltaPlusPlus];
144       fermiEnergy[DeltaPlus] = vDeltaPlus - separationEnergy[DeltaPlus];
145       fermiEnergy[DeltaZero] = vDeltaZero - separationEnergy[DeltaZero];
146       fermiEnergy[DeltaMinus] = vDeltaMinus - separationEnergy[DeltaMinus];
147       
148       fermiEnergy[Lambda] = vLambda - separationEnergy[Lambda];
149       if (fermiEnergy[Lambda] <= 0.)
150          fermiMomentum[Lambda]=0.;
151       else
152          fermiMomentum[Lambda]=std::sqrt(std::pow(fermiEnergy[Lambda]+ml,2.0)-ml*ml);
153 
154       fermiEnergy[SigmaPlus] = vSigmaPlus - separationEnergy[SigmaPlus];
155       fermiEnergy[SigmaZero] = vSigmaZero - separationEnergy[SigmaZero];
156       fermiEnergy[SigmaMinus] = vSigmaMinus - separationEnergy[SigmaMinus];
157    
158       fermiEnergy[antiProton] = vantiProton - separationEnergy[antiProton];
159 
160       INCL_DEBUG("Table of separation energies [MeV] for A=" << theA << ", Z=" << theZ << ":" << '\n'
161             << "  proton:  " << separationEnergy[Proton] << '\n'
162             << "  neutron: " << separationEnergy[Neutron] << '\n'
163             << "  delta++: " << separationEnergy[DeltaPlusPlus] << '\n'
164             << "  delta+:  " << separationEnergy[DeltaPlus] << '\n'
165             << "  delta0:  " << separationEnergy[DeltaZero] << '\n'
166             << "  delta-:  " << separationEnergy[DeltaMinus] << '\n'
167             << "  pi+:     " << separationEnergy[PiPlus] << '\n'
168             << "  pi0:     " << separationEnergy[PiZero] << '\n'
169             << "  pi-:     " << separationEnergy[PiMinus] << '\n'
170             << "  eta:     " << separationEnergy[Eta] << '\n'
171             << "  omega:   " << separationEnergy[Omega] << '\n'
172             << "  etaprime:" << separationEnergy[EtaPrime] << '\n'
173             << "  photon:  " << separationEnergy[Photon] << '\n'
174             << "  lambda:  " << separationEnergy[Lambda] << '\n'
175             << "  sigmaplus:  " << separationEnergy[SigmaPlus] << '\n'
176             << "  sigmazero:  " << separationEnergy[SigmaZero] << '\n'
177             << "  sigmaminus:  " << separationEnergy[SigmaMinus] << '\n'
178             << "  kplus:  " << separationEnergy[KPlus] << '\n'
179             << "  kzero:  " << separationEnergy[KZero] << '\n'
180             << "  kzerobar:  " << separationEnergy[KZeroBar] << '\n'
181             << "  kminus:  " << separationEnergy[KMinus] << '\n'
182             << "  kshort:  " << separationEnergy[KShort] << '\n'
183             << "  klong:  " << separationEnergy[KLong] << '\n'
184             );
185 
186       INCL_DEBUG("Table of Fermi energies [MeV] for A=" << theA << ", Z=" << theZ << ":" << '\n'
187             << "  proton:  " << fermiEnergy[Proton] << '\n'
188             << "  neutron: " << fermiEnergy[Neutron] << '\n'
189             << "  delta++: " << fermiEnergy[DeltaPlusPlus] << '\n'
190             << "  delta+:  " << fermiEnergy[DeltaPlus] << '\n'
191             << "  delta0:  " << fermiEnergy[DeltaZero] << '\n'
192             << "  delta-:  " << fermiEnergy[DeltaMinus] << '\n'
193             << "  lambda:  " << fermiEnergy[Lambda] << '\n'
194             << "  sigma+:  " << fermiEnergy[SigmaPlus] << '\n'
195             << "  sigma0:  " << fermiEnergy[SigmaZero] << '\n'
196             << "  sigma-:  " << fermiEnergy[SigmaMinus] << '\n'
197             );
198 
199       INCL_DEBUG("Table of Fermi momenta [MeV/c] for A=" << theA << ", Z=" << theZ << ":" << '\n'
200             << "  proton:  " << fermiMomentum[Proton] << '\n'
201             << "  neutron: " << fermiMomentum[Neutron] << '\n'
202             );
203     }
204 
205     G4double NuclearPotentialIsospin::computePotentialEnergy(const Particle *particle) const {
206 
207       switch( particle->getType() )
208       {
209         case Proton:
210           return vProton;
211           break;
212         case Neutron:
213           return vNeutron;
214           break;
215 
216         case PiPlus:
217         case PiZero:
218         case PiMinus:
219           return computePionPotentialEnergy(particle);
220           break;
221         
222         case SigmaPlus:
223           return vSigmaPlus;
224           break;
225         case SigmaZero:
226           return vSigmaZero;
227           break;
228         case Lambda:
229           return vLambda;
230           break;
231         case SigmaMinus:
232           return vSigmaMinus;
233           break;
234 
235         case Eta:
236         case Omega:
237         case EtaPrime:
238           return computePionResonancePotentialEnergy(particle);
239           break;
240 
241         case KPlus:
242         case KZero:
243         case KZeroBar:
244         case KMinus:
245         case KShort:
246         case KLong:
247           return computeKaonPotentialEnergy(particle);
248           break;
249 
250         case Photon:
251           return 0.0;
252           break;
253 
254         case antiProton:
255           return vantiProton;
256           break;
257         case antiNeutron:
258           return vantiProton;
259           break;
260         case antiLambda:
261           return 0.0;
262           break;
263         case antiSigmaMinus:
264           return 0.0;
265           break;
266         case antiSigmaPlus:
267           return 0.0;
268           break;
269         case antiSigmaZero:
270           return 0.0;
271           break;
272         case antiXiMinus:
273           return 0.0;
274           break;
275         case antiXiZero:
276           return 0.0;
277           break;
278         case XiMinus:
279           return 0.0;
280           break;
281         case XiZero:
282           return 0.0;
283           break;
284 
285         case DeltaPlusPlus:
286           return vDeltaPlusPlus;
287           break;
288         case DeltaPlus:
289           return vDeltaPlus;
290           break;
291         case DeltaZero:
292           return vDeltaZero;
293           break;
294         case DeltaMinus:
295           return vDeltaMinus;
296           break;
297       case Composite:
298   INCL_ERROR("No potential computed for particle of type Cluster.");
299   return 0.0;
300   break;
301       case UnknownParticle:
302   INCL_ERROR("Trying to compute potential energy for an unknown particle.");
303   return 0.0;
304   break;
305       }
306 
307       INCL_ERROR("There is no potential for this type of particle.");
308       return 0.0;
309     }
310 
311   }
312 }
313 
314