Geant4 Cross Reference |
1 // 1 2 // ******************************************* 3 // * License and Disclaimer 4 // * 5 // * The Geant4 software is copyright of th 6 // * the Geant4 Collaboration. It is provided 7 // * conditions of the Geant4 Software License 8 // * LICENSE and available at http://cern.ch/ 9 // * include a list of copyright holders. 10 // * 11 // * Neither the authors of this software syst 12 // * institutes,nor the agencies providing fin 13 // * work make any representation or warran 14 // * regarding this software system or assum 15 // * use. Please see the license in the file 16 // * for the full disclaimer and the limitatio 17 // * 18 // * This code implementation is the result 19 // * technical work of the GEANT4 collaboratio 20 // * By using, copying, modifying or distri 21 // * any work based on the software) you ag 22 // * use in resulting scientific publicati 23 // * acceptance of all terms of the Geant4 Sof 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 H 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 poten 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::NuclearPotentialI 58 : INuclearPotential(A, Z, aPionPotential 59 { 60 initialize(); 61 } 62 63 // Destructor 64 NuclearPotentialIsospin::~NuclearPotential 65 66 void NuclearPotentialIsospin::initialize() 67 const G4double ZOverA = ((G4double) theZ 68 69 const G4double mp = ParticleTable::getIN 70 const G4double mn = ParticleTable::getIN 71 const G4double ml = ParticleTable::getIN 72 73 const G4double theFermiMomentum = Partic 74 75 fermiMomentum[Proton] = theFermiMomentum 76 const G4double theProtonFermiEnergy = st 77 fermiEnergy[Proton] = theProtonFermiEner 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 98 const G4double tinyMargin = 1E-7; 99 vDeltaPlus = vProton; 100 vDeltaZero = vNeutron; 101 vDeltaPlusPlus = std::max(separationEner 102 vDeltaMinus = std::max(separationEnergyD 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 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 } 204 205 G4double NuclearPotentialIsospin::computeP 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(pa 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; 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 partic 299 return 0.0; 300 break; 301 case UnknownParticle: 302 INCL_ERROR("Trying to compute potential ener 303 return 0.0; 304 break; 305 } 306 307 INCL_ERROR("There is no potential for th 308 return 0.0; 309 } 310 311 } 312 } 313 314