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 9.5)


  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          <<  32 //
                                                   >>  33 // INCL++ revision: v5.0_rc3
 33 //                                                 34 //
 34 #define INCLXX_IN_GEANT4_MODE 1                    35 #define INCLXX_IN_GEANT4_MODE 1
 35                                                    36 
 36 #include "globals.hh"                              37 #include "globals.hh"
 37                                                    38 
 38 /** \file G4INCLNuclearPotentialIsospin.cc         39 /** \file G4INCLNuclearPotentialIsospin.cc
 39  * \brief Isospin-dependent nuclear potential.     40  * \brief Isospin-dependent nuclear potential.
 40  *                                                 41  *
 41  * Provides an isospin-dependent nuclear poten     42  * Provides an isospin-dependent nuclear potential.
 42  *                                                 43  *
 43  * \date 28 February 2011                      <<  44  * Created on: 28 February 2011
 44  * \author Davide Mancusi                      <<  45  *     Author: Davide Mancusi
 45  */                                                46  */
 46                                                    47 
 47 #include "G4INCLNuclearPotentialIsospin.hh"        48 #include "G4INCLNuclearPotentialIsospin.hh"
 48 #include "G4INCLNuclearPotentialConstant.hh"       49 #include "G4INCLNuclearPotentialConstant.hh"
 49 #include "G4INCLParticleTable.hh"                  50 #include "G4INCLParticleTable.hh"
 50 #include "G4INCLGlobals.hh"                        51 #include "G4INCLGlobals.hh"
 51                                                    52 
 52 namespace G4INCL {                                 53 namespace G4INCL {
 53                                                    54 
 54   namespace NuclearPotential {                     55   namespace NuclearPotential {
 55                                                    56 
 56     // Constructors                                57     // Constructors
 57     NuclearPotentialIsospin::NuclearPotentialI <<  58     NuclearPotentialIsospin::NuclearPotentialIsospin(NuclearDensity *density, G4bool pionPotential)
 58       : INuclearPotential(A, Z, aPionPotential <<  59       : INuclearPotential(density, pionPotential)
 59     {                                              60     {
 60       initialize();                                61       initialize();
 61     }                                              62     }
 62                                                    63 
 63     // Destructor                                  64     // Destructor
 64     NuclearPotentialIsospin::~NuclearPotential     65     NuclearPotentialIsospin::~NuclearPotentialIsospin() {}
 65                                                    66 
 66     void NuclearPotentialIsospin::initialize()     67     void NuclearPotentialIsospin::initialize() {
 67       const G4double ZOverA = ((G4double) theZ <<  68       const G4double ZOverA = ((G4double) theDensity->getZ()) / ((G4double) theDensity->getA());
 68                                                    69 
 69       const G4double mp = ParticleTable::getIN <<  70       const G4double mp = ParticleTable::getMass(Proton);
 70       const G4double mn = ParticleTable::getIN <<  71       fermiMomentum[Proton] = Pf * Math::pow13(2.*ZOverA);
 71       const G4double ml = ParticleTable::getIN <<  72       fermiEnergy[Proton] = std::sqrt(fermiMomentum[Proton]*fermiMomentum[Proton] + mp*mp) - mp;
 72                                                <<  73       vProton = fermiEnergy[Proton] + ParticleTable::getSeparationEnergy(Proton);
 73       const G4double theFermiMomentum = Partic <<  74 
 74                                                <<  75       const G4double mn = ParticleTable::getMass(Neutron);
 75       fermiMomentum[Proton] = theFermiMomentum <<  76       fermiMomentum[Neutron] = Pf * Math::pow13(2.*(1.-ZOverA));
 76       const G4double theProtonFermiEnergy = st <<  77       fermiEnergy[Neutron] = std::sqrt(fermiMomentum[Neutron]*fermiMomentum[Neutron] + mn*mn) - mn;
 77       fermiEnergy[Proton] = theProtonFermiEner <<  78       vNeutron = fermiEnergy[Neutron] + ParticleTable::getSeparationEnergy(Neutron);
 78       // Use separation energies from the Part << 
 79       const G4double theProtonSeparationEnergy << 
 80       separationEnergy[Proton] = theProtonSepa << 
 81       vProton = theProtonFermiEnergy + theProt << 
 82                                                << 
 83       fermiMomentum[Neutron] = theFermiMomentu << 
 84       const G4double theNeutronFermiEnergy = s << 
 85       fermiEnergy[Neutron] = theNeutronFermiEn << 
 86       // Use separation energies from the Part << 
 87       const G4double theNeutronSeparationEnerg << 
 88       separationEnergy[Neutron] = theNeutronSe << 
 89       vNeutron = theNeutronFermiEnergy + theNe << 
 90                                                << 
 91       const G4double separationEnergyDeltaPlus << 
 92       separationEnergy[DeltaPlusPlus] = separa << 
 93       separationEnergy[DeltaPlus] = theProtonS << 
 94       separationEnergy[DeltaZero] = theNeutron << 
 95       const G4double separationEnergyDeltaMinu << 
 96       separationEnergy[DeltaMinus] = separatio << 
 97                                                    79 
 98       const G4double tinyMargin = 1E-7;        << 
 99       vDeltaPlus = vProton;                        80       vDeltaPlus = vProton;
100       vDeltaZero = vNeutron;                       81       vDeltaZero = vNeutron;
101       vDeltaPlusPlus = std::max(separationEner <<  82       vDeltaPlusPlus = 2*vDeltaPlus - vDeltaZero;
102       vDeltaMinus = std::max(separationEnergyD <<  83       vDeltaMinus = 2*vDeltaZero - vDeltaPlus;
103                                                <<  84 
104       vSigmaMinus = -16.; // Repulsive potenti <<  85       const G4double Tfpp = vDeltaPlusPlus - vProton + fermiEnergy.find(Proton)->second;
105       vSigmaZero = -16.;  // hypothesis: same  <<  86       const G4double Tfp = fermiEnergy.find(Proton)->second;
106       vSigmaPlus = -16.;                       <<  87       const G4double Tf0 = fermiEnergy.find(Neutron)->second;
107                                                <<  88       const G4double Tfm = vDeltaMinus - vNeutron + fermiEnergy.find(Neutron)->second;
108       vLambda = 30.;                           <<  89       fermiEnergy[DeltaPlusPlus] = Tfpp;
109       vantiProton = 100.;                      <<  90       fermiEnergy[DeltaPlus] = Tfp;
110                                                <<  91       fermiEnergy[DeltaZero] = Tf0;
111       const G4double asy = (theA - 2.*theZ)/th <<  92       fermiEnergy[DeltaMinus] = Tfm;
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                                                << 
119       separationEnergy[PiPlus] = theProtonSepa << 
120       separationEnergy[PiZero] = 0.;           << 
121       separationEnergy[PiMinus] = theNeutronSe << 
122                                                << 
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 << 
144       fermiEnergy[DeltaPlus] = vDeltaPlus - se << 
145       fermiEnergy[DeltaZero] = vDeltaZero - se << 
146       fermiEnergy[DeltaMinus] = vDeltaMinus -  << 
147                                                << 
148       fermiEnergy[Lambda] = vLambda - separati << 
149       if (fermiEnergy[Lambda] <= 0.)           << 
150          fermiMomentum[Lambda]=0.;             << 
151       else                                     << 
152          fermiMomentum[Lambda]=std::sqrt(std:: << 
153                                                << 
154       fermiEnergy[SigmaPlus] = vSigmaPlus - se << 
155       fermiEnergy[SigmaZero] = vSigmaZero - se << 
156       fermiEnergy[SigmaMinus] = vSigmaMinus -  << 
157                                                << 
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             );                                 << 
185                                                << 
186       INCL_DEBUG("Table of Fermi energies [MeV << 
187             << "  proton:  " << fermiEnergy[Pr << 
188             << "  neutron: " << fermiEnergy[Ne << 
189             << "  delta++: " << fermiEnergy[De << 
190             << "  delta+:  " << fermiEnergy[De << 
191             << "  delta0:  " << fermiEnergy[De << 
192             << "  delta-:  " << fermiEnergy[De << 
193             << "  lambda:  " << fermiEnergy[La << 
194             << "  sigma+:  " << fermiEnergy[Si << 
195             << "  sigma0:  " << fermiEnergy[Si << 
196             << "  sigma-:  " << fermiEnergy[Si << 
197             );                                 << 
198                                                << 
199       INCL_DEBUG("Table of Fermi momenta [MeV/ << 
200             << "  proton:  " << fermiMomentum[ << 
201             << "  neutron: " << fermiMomentum[ << 
202             );                                 << 
203     }                                              93     }
204                                                    94 
205     G4double NuclearPotentialIsospin::computeP     95     G4double NuclearPotentialIsospin::computePotentialEnergy(const Particle *particle) const {
206                                                    96 
207       switch( particle->getType() )                97       switch( particle->getType() )
208       {                                            98       {
209         case Proton:                               99         case Proton:
210           return vProton;                         100           return vProton;
211           break;                                  101           break;
212         case Neutron:                             102         case Neutron:
213           return vNeutron;                        103           return vNeutron;
214           break;                                  104           break;
215                                                   105 
216         case PiPlus:                              106         case PiPlus:
217         case PiZero:                              107         case PiZero:
218         case PiMinus:                             108         case PiMinus:
219           return computePionPotentialEnergy(pa    109           return computePionPotentialEnergy(particle);
220           break;                                  110           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;                               << 
284                                                   111 
285         case DeltaPlusPlus:                       112         case DeltaPlusPlus:
286           return vDeltaPlusPlus;                  113           return vDeltaPlusPlus;
287           break;                                  114           break;
288         case DeltaPlus:                           115         case DeltaPlus:
289           return vDeltaPlus;                      116           return vDeltaPlus;
290           break;                                  117           break;
291         case DeltaZero:                           118         case DeltaZero:
292           return vDeltaZero;                      119           return vDeltaZero;
293           break;                                  120           break;
294         case DeltaMinus:                          121         case DeltaMinus:
295           return vDeltaMinus;                     122           return vDeltaMinus;
296           break;                                  123           break;
297       case Composite:                             124       case Composite:
298   INCL_ERROR("No potential computed for partic << 125   ERROR("No potential computed for particle of type Cluster.");
299   return 0.0;                                     126   return 0.0;
300   break;                                          127   break;
301       case UnknownParticle:                       128       case UnknownParticle:
302   INCL_ERROR("Trying to compute potential ener << 129   ERROR("Trying to compute potential energy for an unknown particle.");
303   return 0.0;                                     130   return 0.0;
304   break;                                          131   break;
305       }                                           132       }
306                                                   133 
307       INCL_ERROR("There is no potential for th << 134       ERROR("There is no potential for this type of particle.");
308       return 0.0;                                 135       return 0.0;
309     }                                             136     }
310                                                   137 
311   }                                               138   }
312 }                                                 139 }
313                                                   140 
314                                                   141