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 ]

Diff markup

Differences between /processes/hadronic/models/inclxx/incl_physics/src/G4INCLNuclearPotentialIsospin.cc (Version 11.3.0) and /processes/hadronic/models/inclxx/incl_physics/src/G4INCLNuclearPotentialIsospin.cc (Version 10.0.p4)


  1 //                                                  1 //
  2 // *******************************************      2 // ********************************************************************
  3 // * License and Disclaimer                         3 // * License and Disclaimer                                           *
  4 // *                                                4 // *                                                                  *
  5 // * The  Geant4 software  is  copyright of th      5 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
  6 // * the Geant4 Collaboration.  It is provided      6 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
  7 // * conditions of the Geant4 Software License      7 // * conditions of the Geant4 Software License,  included in the file *
  8 // * LICENSE and available at  http://cern.ch/      8 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
  9 // * include a list of copyright holders.           9 // * include a list of copyright holders.                             *
 10 // *                                               10 // *                                                                  *
 11 // * Neither the authors of this software syst     11 // * Neither the authors of this software system, nor their employing *
 12 // * institutes,nor the agencies providing fin     12 // * institutes,nor the agencies providing financial support for this *
 13 // * work  make  any representation or  warran     13 // * work  make  any representation or  warranty, express or implied, *
 14 // * regarding  this  software system or assum     14 // * regarding  this  software system or assume any liability for its *
 15 // * use.  Please see the license in the file      15 // * use.  Please see the license in the file  LICENSE  and URL above *
 16 // * for the full disclaimer and the limitatio     16 // * for the full disclaimer and the limitation of liability.         *
 17 // *                                               17 // *                                                                  *
 18 // * This  code  implementation is the result      18 // * This  code  implementation is the result of  the  scientific and *
 19 // * technical work of the GEANT4 collaboratio     19 // * technical work of the GEANT4 collaboration.                      *
 20 // * By using,  copying,  modifying or  distri     20 // * By using,  copying,  modifying or  distributing the software (or *
 21 // * any work based  on the software)  you  ag     21 // * any work based  on the software)  you  agree  to acknowledge its *
 22 // * use  in  resulting  scientific  publicati     22 // * use  in  resulting  scientific  publications,  and indicate your *
 23 // * acceptance of all terms of the Geant4 Sof     23 // * acceptance of all terms of the Geant4 Software license.          *
 24 // *******************************************     24 // ********************************************************************
 25 //                                                 25 //
 26 // INCL++ intra-nuclear cascade model              26 // INCL++ intra-nuclear cascade model
 27 // Alain Boudard, CEA-Saclay, France           <<  27 // Pekka Kaitaniemi, CEA and Helsinki Institute of Physics
 28 // Joseph Cugnon, University of Liege, Belgium <<  28 // Davide Mancusi, CEA
 29 // Jean-Christophe David, CEA-Saclay, France   <<  29 // Alain Boudard, CEA
 30 // Pekka Kaitaniemi, CEA-Saclay, France, and H <<  30 // Sylvie Leray, CEA
 31 // Sylvie Leray, CEA-Saclay, France            <<  31 // Joseph Cugnon, University of Liege
 32 // Davide Mancusi, CEA-Saclay, France          << 
 33 //                                                 32 //
 34 #define INCLXX_IN_GEANT4_MODE 1                    33 #define INCLXX_IN_GEANT4_MODE 1
 35                                                    34 
 36 #include "globals.hh"                              35 #include "globals.hh"
 37                                                    36 
 38 /** \file G4INCLNuclearPotentialIsospin.cc         37 /** \file G4INCLNuclearPotentialIsospin.cc
 39  * \brief Isospin-dependent nuclear potential.     38  * \brief Isospin-dependent nuclear potential.
 40  *                                                 39  *
 41  * Provides an isospin-dependent nuclear poten     40  * Provides an isospin-dependent nuclear potential.
 42  *                                                 41  *
 43  * \date 28 February 2011                          42  * \date 28 February 2011
 44  * \author Davide Mancusi                          43  * \author Davide Mancusi
 45  */                                                44  */
 46                                                    45 
 47 #include "G4INCLNuclearPotentialIsospin.hh"        46 #include "G4INCLNuclearPotentialIsospin.hh"
 48 #include "G4INCLNuclearPotentialConstant.hh"       47 #include "G4INCLNuclearPotentialConstant.hh"
 49 #include "G4INCLParticleTable.hh"                  48 #include "G4INCLParticleTable.hh"
 50 #include "G4INCLGlobals.hh"                        49 #include "G4INCLGlobals.hh"
 51                                                    50 
 52 namespace G4INCL {                                 51 namespace G4INCL {
 53                                                    52 
 54   namespace NuclearPotential {                     53   namespace NuclearPotential {
 55                                                    54 
 56     // Constructors                                55     // Constructors
 57     NuclearPotentialIsospin::NuclearPotentialI     56     NuclearPotentialIsospin::NuclearPotentialIsospin(const G4int A, const G4int Z, const G4bool aPionPotential)
 58       : INuclearPotential(A, Z, aPionPotential     57       : INuclearPotential(A, Z, aPionPotential)
 59     {                                              58     {
 60       initialize();                                59       initialize();
 61     }                                              60     }
 62                                                    61 
 63     // Destructor                                  62     // Destructor
 64     NuclearPotentialIsospin::~NuclearPotential     63     NuclearPotentialIsospin::~NuclearPotentialIsospin() {}
 65                                                    64 
 66     void NuclearPotentialIsospin::initialize()     65     void NuclearPotentialIsospin::initialize() {
 67       const G4double ZOverA = ((G4double) theZ     66       const G4double ZOverA = ((G4double) theZ) / ((G4double) theA);
 68                                                    67 
 69       const G4double mp = ParticleTable::getIN     68       const G4double mp = ParticleTable::getINCLMass(Proton);
 70       const G4double mn = ParticleTable::getIN     69       const G4double mn = ParticleTable::getINCLMass(Neutron);
 71       const G4double ml = ParticleTable::getIN << 
 72                                                    70 
 73       const G4double theFermiMomentum = Partic     71       const G4double theFermiMomentum = ParticleTable::getFermiMomentum(theA,theZ);
 74                                                    72 
 75       fermiMomentum[Proton] = theFermiMomentum     73       fermiMomentum[Proton] = theFermiMomentum * Math::pow13(2.*ZOverA);
 76       const G4double theProtonFermiEnergy = st     74       const G4double theProtonFermiEnergy = std::sqrt(fermiMomentum[Proton]*fermiMomentum[Proton] + mp*mp) - mp;
 77       fermiEnergy[Proton] = theProtonFermiEner     75       fermiEnergy[Proton] = theProtonFermiEnergy;
 78       // Use separation energies from the Part     76       // Use separation energies from the ParticleTable
 79       const G4double theProtonSeparationEnergy     77       const G4double theProtonSeparationEnergy = ParticleTable::getSeparationEnergy(Proton,theA,theZ);
 80       separationEnergy[Proton] = theProtonSepa     78       separationEnergy[Proton] = theProtonSeparationEnergy;
 81       vProton = theProtonFermiEnergy + theProt     79       vProton = theProtonFermiEnergy + theProtonSeparationEnergy;
 82                                                    80 
 83       fermiMomentum[Neutron] = theFermiMomentu     81       fermiMomentum[Neutron] = theFermiMomentum * Math::pow13(2.*(1.-ZOverA));
 84       const G4double theNeutronFermiEnergy = s     82       const G4double theNeutronFermiEnergy = std::sqrt(fermiMomentum[Neutron]*fermiMomentum[Neutron] + mn*mn) - mn;
 85       fermiEnergy[Neutron] = theNeutronFermiEn     83       fermiEnergy[Neutron] = theNeutronFermiEnergy;
 86       // Use separation energies from the Part     84       // Use separation energies from the ParticleTable
 87       const G4double theNeutronSeparationEnerg     85       const G4double theNeutronSeparationEnergy = ParticleTable::getSeparationEnergy(Neutron,theA,theZ);
 88       separationEnergy[Neutron] = theNeutronSe     86       separationEnergy[Neutron] = theNeutronSeparationEnergy;
 89       vNeutron = theNeutronFermiEnergy + theNe     87       vNeutron = theNeutronFermiEnergy + theNeutronSeparationEnergy;
 90                                                    88 
 91       const G4double separationEnergyDeltaPlus     89       const G4double separationEnergyDeltaPlusPlus = 2.*theProtonSeparationEnergy - theNeutronSeparationEnergy;
 92       separationEnergy[DeltaPlusPlus] = separa     90       separationEnergy[DeltaPlusPlus] = separationEnergyDeltaPlusPlus;
 93       separationEnergy[DeltaPlus] = theProtonS     91       separationEnergy[DeltaPlus] = theProtonSeparationEnergy;
 94       separationEnergy[DeltaZero] = theNeutron     92       separationEnergy[DeltaZero] = theNeutronSeparationEnergy;
 95       const G4double separationEnergyDeltaMinu     93       const G4double separationEnergyDeltaMinus = 2.*theNeutronSeparationEnergy - theProtonSeparationEnergy;
 96       separationEnergy[DeltaMinus] = separatio     94       separationEnergy[DeltaMinus] = separationEnergyDeltaMinus;
 97                                                    95 
 98       const G4double tinyMargin = 1E-7;            96       const G4double tinyMargin = 1E-7;
 99       vDeltaPlus = vProton;                        97       vDeltaPlus = vProton;
100       vDeltaZero = vNeutron;                       98       vDeltaZero = vNeutron;
101       vDeltaPlusPlus = std::max(separationEner     99       vDeltaPlusPlus = std::max(separationEnergyDeltaPlusPlus + tinyMargin, 2.*vDeltaPlus - vDeltaZero);
102       vDeltaMinus = std::max(separationEnergyD    100       vDeltaMinus = std::max(separationEnergyDeltaMinus + tinyMargin, 2.*vDeltaZero - vDeltaPlus);
103                                                << 
104       vSigmaMinus = -16.; // Repulsive potenti << 
105       vSigmaZero = -16.;  // hypothesis: same  << 
106       vSigmaPlus = -16.;                       << 
107                                                << 
108       vLambda = 30.;                           << 
109       vantiProton = 100.;                      << 
110                                                << 
111       const G4double asy = (theA - 2.*theZ)/th << 
112       // Jose Luis Rodriguez-Sanchez et al., R << 
113       if      (asy > 0.236) vLambda = 40.91;   << 
114       else if (asy > 0.133) vLambda = 56.549 - << 
115                                                << 
116       const G4double theLambdaSeparationEnergy << 
117       const G4double theantiProtonSeparationEn << 
118                                                   101 
119       separationEnergy[PiPlus] = theProtonSepa    102       separationEnergy[PiPlus] = theProtonSeparationEnergy - theNeutronSeparationEnergy;
120       separationEnergy[PiZero] = 0.;              103       separationEnergy[PiZero] = 0.;
121       separationEnergy[PiMinus] = theNeutronSe    104       separationEnergy[PiMinus] = theNeutronSeparationEnergy - theProtonSeparationEnergy;
122                                                   105 
123       separationEnergy[Eta]      = 0.;         << 
124       separationEnergy[Omega]    = 0.;         << 
125       separationEnergy[EtaPrime] = 0.;         << 
126       separationEnergy[Photon]   = 0.;         << 
127                                                << 
128       separationEnergy[Lambda]    = theLambdaS << 
129       separationEnergy[SigmaPlus] = theProtonS << 
130       separationEnergy[SigmaZero] = theLambdaS << 
131       separationEnergy[SigmaMinus]  = theNeutr << 
132                                                << 
133       separationEnergy[KPlus]   = theProtonSep << 
134       separationEnergy[KZero]   = (theNeutronS << 
135       separationEnergy[KZeroBar]  = (theLambda << 
136       separationEnergy[KMinus]    = 2.*theNeut << 
137                                                << 
138       separationEnergy[KShort]    = (theNeutro << 
139       separationEnergy[KLong]   = (theNeutronS << 
140                                                << 
141       separationEnergy[antiProton]    = theant << 
142                                                << 
143       fermiEnergy[DeltaPlusPlus] = vDeltaPlusP    106       fermiEnergy[DeltaPlusPlus] = vDeltaPlusPlus - separationEnergy[DeltaPlusPlus];
144       fermiEnergy[DeltaPlus] = vDeltaPlus - se    107       fermiEnergy[DeltaPlus] = vDeltaPlus - separationEnergy[DeltaPlus];
145       fermiEnergy[DeltaZero] = vDeltaZero - se    108       fermiEnergy[DeltaZero] = vDeltaZero - separationEnergy[DeltaZero];
146       fermiEnergy[DeltaMinus] = vDeltaMinus -     109       fermiEnergy[DeltaMinus] = vDeltaMinus - separationEnergy[DeltaMinus];
147                                                << 110 
148       fermiEnergy[Lambda] = vLambda - separati << 111       INCL_DEBUG("Table of separation energies [MeV] for A=" << theA << ", Z=" << theZ << ":" << std::endl
149       if (fermiEnergy[Lambda] <= 0.)           << 112             << "  proton:  " << separationEnergy[Proton] << std::endl
150          fermiMomentum[Lambda]=0.;             << 113             << "  neutron: " << separationEnergy[Neutron] << std::endl
151       else                                     << 114             << "  delta++: " << separationEnergy[DeltaPlusPlus] << std::endl
152          fermiMomentum[Lambda]=std::sqrt(std:: << 115             << "  delta+:  " << separationEnergy[DeltaPlus] << std::endl
153                                                << 116             << "  delta0:  " << separationEnergy[DeltaZero] << std::endl
154       fermiEnergy[SigmaPlus] = vSigmaPlus - se << 117             << "  delta-:  " << separationEnergy[DeltaMinus] << std::endl
155       fermiEnergy[SigmaZero] = vSigmaZero - se << 118             << "  pi+:     " << separationEnergy[PiPlus] << std::endl
156       fermiEnergy[SigmaMinus] = vSigmaMinus -  << 119             << "  pi0:     " << separationEnergy[PiZero] << std::endl
157                                                << 120             << "  pi-:     " << separationEnergy[PiMinus] << std::endl
158       fermiEnergy[antiProton] = vantiProton -  << 
159                                                << 
160       INCL_DEBUG("Table of separation energies << 
161             << "  proton:  " << separationEner << 
162             << "  neutron: " << separationEner << 
163             << "  delta++: " << separationEner << 
164             << "  delta+:  " << separationEner << 
165             << "  delta0:  " << separationEner << 
166             << "  delta-:  " << separationEner << 
167             << "  pi+:     " << separationEner << 
168             << "  pi0:     " << separationEner << 
169             << "  pi-:     " << separationEner << 
170             << "  eta:     " << separationEner << 
171             << "  omega:   " << separationEner << 
172             << "  etaprime:" << separationEner << 
173             << "  photon:  " << separationEner << 
174             << "  lambda:  " << separationEner << 
175             << "  sigmaplus:  " << separationE << 
176             << "  sigmazero:  " << separationE << 
177             << "  sigmaminus:  " << separation << 
178             << "  kplus:  " << separationEnerg << 
179             << "  kzero:  " << separationEnerg << 
180             << "  kzerobar:  " << separationEn << 
181             << "  kminus:  " << separationEner << 
182             << "  kshort:  " << separationEner << 
183             << "  klong:  " << separationEnerg << 
184             );                                    121             );
185                                                   122 
186       INCL_DEBUG("Table of Fermi energies [MeV << 123       INCL_DEBUG("Table of Fermi energies [MeV] for A=" << theA << ", Z=" << theZ << ":" << std::endl
187             << "  proton:  " << fermiEnergy[Pr << 124             << "  proton:  " << fermiEnergy[Proton] << std::endl
188             << "  neutron: " << fermiEnergy[Ne << 125             << "  neutron: " << fermiEnergy[Neutron] << std::endl
189             << "  delta++: " << fermiEnergy[De << 126             << "  delta++: " << fermiEnergy[DeltaPlusPlus] << std::endl
190             << "  delta+:  " << fermiEnergy[De << 127             << "  delta+:  " << fermiEnergy[DeltaPlus] << std::endl
191             << "  delta0:  " << fermiEnergy[De << 128             << "  delta0:  " << fermiEnergy[DeltaZero] << std::endl
192             << "  delta-:  " << fermiEnergy[De << 129             << "  delta-:  " << fermiEnergy[DeltaMinus] << std::endl
193             << "  lambda:  " << fermiEnergy[La << 
194             << "  sigma+:  " << fermiEnergy[Si << 
195             << "  sigma0:  " << fermiEnergy[Si << 
196             << "  sigma-:  " << fermiEnergy[Si << 
197             );                                    130             );
198                                                   131 
199       INCL_DEBUG("Table of Fermi momenta [MeV/ << 132       INCL_DEBUG("Table of Fermi momenta [MeV/c] for A=" << theA << ", Z=" << theZ << ":" << std::endl
200             << "  proton:  " << fermiMomentum[ << 133             << "  proton:  " << fermiMomentum[Proton] << std::endl
201             << "  neutron: " << fermiMomentum[ << 134             << "  neutron: " << fermiMomentum[Neutron] << std::endl
202             );                                    135             );
203     }                                             136     }
204                                                   137 
205     G4double NuclearPotentialIsospin::computeP    138     G4double NuclearPotentialIsospin::computePotentialEnergy(const Particle *particle) const {
206                                                   139 
207       switch( particle->getType() )               140       switch( particle->getType() )
208       {                                           141       {
209         case Proton:                              142         case Proton:
210           return vProton;                         143           return vProton;
211           break;                                  144           break;
212         case Neutron:                             145         case Neutron:
213           return vNeutron;                        146           return vNeutron;
214           break;                                  147           break;
215                                                   148 
216         case PiPlus:                              149         case PiPlus:
217         case PiZero:                              150         case PiZero:
218         case PiMinus:                             151         case PiMinus:
219           return computePionPotentialEnergy(pa    152           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 computePionResonancePotential << 
239           break;                               << 
240                                                << 
241         case KPlus:                            << 
242         case KZero:                            << 
243         case KZeroBar:                         << 
244         case KMinus:                           << 
245         case KShort:                           << 
246         case KLong:                            << 
247           return computeKaonPotentialEnergy(pa << 
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;                                  153           break;
284                                                   154 
285         case DeltaPlusPlus:                       155         case DeltaPlusPlus:
286           return vDeltaPlusPlus;                  156           return vDeltaPlusPlus;
287           break;                                  157           break;
288         case DeltaPlus:                           158         case DeltaPlus:
289           return vDeltaPlus;                      159           return vDeltaPlus;
290           break;                                  160           break;
291         case DeltaZero:                           161         case DeltaZero:
292           return vDeltaZero;                      162           return vDeltaZero;
293           break;                                  163           break;
294         case DeltaMinus:                          164         case DeltaMinus:
295           return vDeltaMinus;                     165           return vDeltaMinus;
296           break;                                  166           break;
297       case Composite:                             167       case Composite:
298   INCL_ERROR("No potential computed for partic    168   INCL_ERROR("No potential computed for particle of type Cluster.");
299   return 0.0;                                     169   return 0.0;
300   break;                                          170   break;
301       case UnknownParticle:                       171       case UnknownParticle:
302   INCL_ERROR("Trying to compute potential ener    172   INCL_ERROR("Trying to compute potential energy for an unknown particle.");
303   return 0.0;                                     173   return 0.0;
304   break;                                          174   break;
305       }                                           175       }
306                                                   176 
307       INCL_ERROR("There is no potential for th    177       INCL_ERROR("There is no potential for this type of particle.");
308       return 0.0;                                 178       return 0.0;
309     }                                             179     }
310                                                   180 
311   }                                               181   }
312 }                                                 182 }
313                                                   183 
314                                                   184