include/Geant4/G4INCLINuclearPotential.hh

0001 //
0002 // ********************************************************************
0003 // * License and Disclaimer                                           *
0004 // *                                                                  *
0005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
0006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
0007 // * conditions of the Geant4 Software License,  included in the file *
0008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
0009 // * include a list of copyright holders.                             *
0010 // *                                                                  *
0011 // * Neither the authors of this software system, nor their employing *
0012 // * institutes,nor the agencies providing financial support for this *
0013 // * work  make  any representation or  warranty, express or implied, *
0014 // * regarding  this  software system or assume any liability for its *
0015 // * use.  Please see the license in the file  LICENSE  and URL above *
0016 // * for the full disclaimer and the limitation of liability.         *
0017 // *                                                                  *
0018 // * This  code  implementation is the result of  the  scientific and *
0019 // * technical work of the GEANT4 collaboration.                      *
0020 // * By using,  copying,  modifying or  distributing the software (or *
0021 // * any work based  on the software)  you  agree  to acknowledge its *
0022 // * use  in  resulting  scientific  publications,  and indicate your *
0023 // * acceptance of all terms of the Geant4 Software license.          *
0024 // ********************************************************************
0025 //
0026 // INCL++ intra-nuclear cascade model
0027 // Alain Boudard, CEA-Saclay, France
0028 // Joseph Cugnon, University of Liege, Belgium
0029 // Jean-Christophe David, CEA-Saclay, France
0030 // Pekka Kaitaniemi, CEA-Saclay, France, and Helsinki Institute of Physics, Finland
0031 // Sylvie Leray, CEA-Saclay, France
0032 // Davide Mancusi, CEA-Saclay, France
0033 //
0034 #define INCLXX_IN_GEANT4_MODE 1
0035
0036 #include "globals.hh"
0037
0038 /** \file G4INCLINuclearPotential.hh
0039  * \brief Abstract interface to the nuclear potential.
0040  *
0041  * NuclearPotential-like classes should provide access to the value of the
0042  * potential of a particle in a particular context. For example, an instance of
0043  * a NuclearPotential class should be associated to every nucleus.
0044  *
0045  * \date 17 January 2011
0046  * \author Davide Mancusi
0047  */
0048
0049 #ifndef G4INCLINUCLEARPOTENTIAL_HH
0050 #define G4INCLINUCLEARPOTENTIAL_HH 1
0051
0052 #include "G4INCLParticle.hh"
0053 #include "G4INCLRandom.hh"
0054 #include "G4INCLDeuteronDensity.hh"
0055 #include <map>
0056 // #include <cassert>
0057
0058 namespace G4INCL {
0059
0060   namespace NuclearPotential {
0061     class INuclearPotential {
0062       public:
0063         INuclearPotential(const G4int A, const G4int Z, const G4bool pionPot) :
0064           theA(A),
0065           theZ(Z),
0066           pionPotential(pionPot)
0067         {
0068           if(pionPotential) {
0069             const G4double ZOverA = ((G4double) theZ) / ((G4double) theA);
0070             // As in INCL4.6, use the r0*A^(1/3) formula to estimate vc
0071             const G4double r = 1.12*Math::pow13((G4double)theA);
0072
0073             const G4double xsi = 1. - 2.*ZOverA;
0074             const G4double vc = 1.25*PhysicalConstants::eSquared*theZ/r;
0075             vPiPlus = vPionDefault + 71.*xsi - vc;
0076             vPiZero = vPionDefault;
0077             vPiMinus = vPionDefault - 71.*xsi + vc;
0078             vKPlus = vKPlusDefault;
0079             vKZero = vKPlusDefault + 10.; // Hypothesis to be check
0080             vKMinus = vKMinusDefault;
0081             vKZeroBar = vKMinusDefault - 10.; // Hypothesis to be check
0082           } else {
0083             vPiPlus = 0.0;
0084             vPiZero = 0.0;
0085             vPiMinus = 0.0;
0086             vKPlus = 0.0;
0087             vKZero = 0.0;
0088             vKMinus = 0.0;
0089             vKZeroBar = 0.0;
0090           }
0091         }
0092
0093         virtual ~INuclearPotential() {}
0094
0095         /// \brief Do we have a pion potential?
0096         G4bool hasPionPotential() const { return pionPotential; }
0097
0098         virtual G4double computePotentialEnergy(const Particle * const p) const = 0;
0099
0100         /** \brief Return the Fermi energy for a particle.
0101          *
0102          * \param p pointer to a Particle
0103          * \return Fermi energy for that particle type
0104          **/
0105         inline G4double getFermiEnergy(const Particle * const p) const {
0106           std::map<ParticleType, G4double>::const_iterator i = fermiEnergy.find(p->getType());
0107 // assert(i!=fermiEnergy.end());
0108           return i->second;
0109         }
0110
0111         /** \brief Return the Fermi energy for a particle type.
0112          *
0113          * \param t particle type
0114          * \return Fermi energy for that particle type
0115          **/
0116         inline G4double getFermiEnergy(const ParticleType t) const {
0117           std::map<ParticleType, G4double>::const_iterator i = fermiEnergy.find(t);
0118 // assert(i!=fermiEnergy.end());
0119           return i->second;
0120         }
0121
0122         /** \brief Return the separation energy for a particle.
0123          *
0124          * \param p pointer to a Particle
0125          * \return separation energy for that particle type
0126          **/
0127         inline G4double getSeparationEnergy(const Particle * const p) const {
0128           std::map<ParticleType, G4double>::const_iterator i = separationEnergy.find(p->getType());
0129 // assert(i!=separationEnergy.end());
0130           return i->second;
0131         }
0132
0133         /** \brief Return the separation energy for a particle type.
0134          *
0135          * \param t particle type
0136          * \return separation energy for that particle type
0137          **/
0138         inline G4double getSeparationEnergy(const ParticleType t) const {
0139           std::map<ParticleType, G4double>::const_iterator i = separationEnergy.find(t);
0140 // assert(i!=separationEnergy.end());
0141           return i->second;
0142         }
0143
0144         /** \brief Return the Fermi momentum for a particle.
0145          *
0146          * \param p pointer to a Particle
0147          * \return Fermi momentum for that particle type
0148          **/
0149         inline G4double getFermiMomentum(const Particle * const p) const {
0150           if(p->isDelta()) {
0151             const G4double Tf = getFermiEnergy(p), mass = p->getMass();
0152             return std::sqrt(Tf*(Tf+2.*mass));
0153           } else {
0154             std::map<ParticleType, G4double>::const_iterator i = fermiMomentum.find(p->getType());
0155 // assert(i!=fermiMomentum.end());
0156             return i->second;
0157           }
0158         }
0159
0160         /** \brief Return the Fermi momentum for a particle type.
0161          *
0162          * \param t particle type
0163          * \return Fermi momentum for that particle type
0164          **/
0165         inline G4double getFermiMomentum(const ParticleType t) const {
0166 // assert(t!=DeltaPlusPlus && t!=DeltaPlus && t!=DeltaZero && t!=DeltaMinus);
0167           std::map<ParticleType, G4double>::const_iterator i = fermiMomentum.find(t);
0168           return i->second;
0169         }
0170
0171       protected:
0172         /// \brief Compute the potential energy for the given pion.
0173         G4double computePionPotentialEnergy(const Particle * const p) const {
0174 // assert(p->getType()==PiPlus || p->getType()==PiZero || p->getType()==PiMinus);
0175           if(pionPotential && !p->isOutOfWell()) {
0176             switch( p->getType() ) {
0177               case PiPlus:
0178                 return vPiPlus;
0179                 break;
0180               case PiZero:
0181                 return vPiZero;
0182                 break;
0183               case PiMinus:
0184                 return vPiMinus;
0185                 break;
0186               default: // Pion potential is defined and non-zero only for pions
0187                 return 0.0;
0188                 break;
0189             }
0190           }
0191           else
0192             return 0.0;
0193         }
0194
0195       protected:
0196         /// \brief Compute the potential energy for the given kaon.
0197         G4double computeKaonPotentialEnergy(const Particle * const p) const {
0198 // assert(p->getType()==KPlus || p->getType()==KZero || p->getType()==KZeroBar || p->getType()==KMinus|| p->getType()==KShort|| p->getType()==KLong);
0199           if(pionPotential && !p->isOutOfWell()) { // if pionPotental false -> kaonPotential false
0200             switch( p->getType() ) {
0201               case KPlus:
0202                 return vKPlus;
0203                 break;
0204               case KZero:
0205                 return vKZero;
0206                 break;
0207               case KZeroBar:
0208                 return vKZeroBar;
0209                 break;
0210               case KShort:
0211               case KLong:
0212                  return 0.0; // Should never be in the nucleus
0213                  break;
0214                case KMinus:
0215                  return vKMinus;
0216                  break;
0217                default:
0218                  return 0.0;
0219                  break;
0220                }
0221             }
0222         else
0223           return 0.0;
0224         }
0225
0226       protected:
0227         /// \brief Compute the potential energy for the given pion resonances (Eta, Omega and EtaPrime and Gamma also).
0228         G4double computePionResonancePotentialEnergy(const Particle * const p) const {
0229 // assert(p->getType()==Eta || p->getType()==Omega || p->getType()==EtaPrime || p->getType()==Photon);
0230           if(pionPotential && !p->isOutOfWell()) {
0231             switch( p->getType() ) {
0232               case Eta:
0233 //jcd             return vPiZero;
0234 //jcd              return vPiZero*1.5;
0235                 return 0.0; // (JCD: seems to give better results)
0236                 break;
0237               case Omega:
0238                 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)
0239                 break;
0240               case EtaPrime:
0241                 return 37.0; // V. Metag in Hyperfine Interact (2015) 234:25-31
0242                 break;
0243               case Photon:
0244                 return 0.0;
0245                 break;
0246               default:
0247                 return 0.0;
0248                 break;
0249             }
0250           }
0251           else
0252             return 0.0;
0253         }
0254
0255       protected:
0256         /// \brief The mass number of the nucleus
0257         const G4int theA;
0258         /// \brief The charge number of the nucleus
0259         const G4int theZ;
0260       private:
0261         const G4bool pionPotential;
0262         G4double vPiPlus, vPiZero, vPiMinus;
0263         static const G4double vPionDefault;
0264         G4double vKPlus, vKZero, vKZeroBar, vKMinus;
0265         static const G4double vKPlusDefault;
0266         static const G4double vKMinusDefault;
0267       protected:
0268         /* \brief map of Fermi energies per particle type */
0269         std::map<ParticleType,G4double> fermiEnergy;
0270         /* \brief map of Fermi momenta per particle type */
0271         std::map<ParticleType,G4double> fermiMomentum;
0272         /* \brief map of separation energies per particle type */
0273         std::map<ParticleType,G4double> separationEnergy;
0274
0275     };
0276
0277
0278
0279     /** \brief Create an INuclearPotential object
0280      *
0281      * This is the method that should be used to instantiate objects derived
0282      * from INuclearPotential. It uses a caching mechanism to minimise
0283      * thrashing and speed up the code.
0284      *
0285      * \param type the type of the potential to be created
0286      * \param theA mass number of the nucleus
0287      * \param theZ charge number of the nucleus
0288      * \param pionPotential whether pions should also feel the potential
0289      * \return a pointer to the nuclear potential
0290      */
0291     INuclearPotential const *createPotential(const PotentialType type, const G4int theA, const G4int theZ, const G4bool pionPotential);
0292
0293     /// \brief Clear the INuclearPotential cache
0294     void clearCache();
0295
0296   }
0297
0298 }
0299
0300 #endif /* G4INCLINUCLEARPOTENTIAL_HH_ */