Geant4 Cross Reference |
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 // 26 // 27 // 27 // 28 // original by H.P. Wellisch << 28 // original by H.P. Wellisch 29 // modified by J.L. Chuma, TRIUMF, 19-Nov-1996 << 29 // modified by J.L. Chuma, TRIUMF, 19-Nov-1996 30 // last modified: 27-Mar-1997 << 30 // last modified: 27-Mar-1997 31 // J.P.Wellisch: 23-Apr-97: minor simplificati << 31 // J.P.Wellisch: 23-Apr-97: minor simplifications 32 // modified by J.L.Chuma 24-Jul-97 to set the << 32 // modified by J.L.Chuma 24-Jul-97 to set the total momentum in Cinema and 33 // Evaporatio << 33 // EvaporationEffects 34 // modified by J.L.Chuma 21-Oct-97 put std::a << 34 // modified by J.L.Chuma 21-Oct-97 put std::abs() around the totalE^2-mass^2 35 // in calcula << 35 // in calculation of total momentum in 36 // Cinema and << 36 // Cinema and EvaporationEffects 37 // Chr. Volcker, 10-Nov-1997: new methods and << 37 // Chr. Volcker, 10-Nov-1997: new methods and class variables. 38 // HPW added utilities for low energy neutron << 38 // HPW added utilities for low energy neutron transport. (12.04.1998) 39 // M.G. Pia, 2 Oct 1998: modified GetFermiMome << 39 // M.G. Pia, 2 Oct 1998: modified GetFermiMomentum to avoid memory leaks 40 // G.Folger, spring 2010: add integer A/Z int << 40 // G.Folger, spring 2010: add integer A/Z interface 41 // A. Ribon, summer 2015: migrated to G4Exp a << 42 // A. Ribon, autumn 2021: extended to hypernu << 43 41 44 #include "G4Nucleus.hh" 42 #include "G4Nucleus.hh" 45 #include "G4NucleiProperties.hh" 43 #include "G4NucleiProperties.hh" 46 #include "G4PhysicalConstants.hh" << 47 #include "G4SystemOfUnits.hh" << 48 #include "Randomize.hh" 44 #include "Randomize.hh" 49 #include "G4HadronicException.hh" 45 #include "G4HadronicException.hh" 50 #include "G4Exp.hh" << 46 51 #include "G4Log.hh" << 52 #include "G4HyperNucleiProperties.hh" << 53 #include "G4HadronicParameters.hh" << 54 << 55 << 56 G4Nucleus::G4Nucleus() 47 G4Nucleus::G4Nucleus() 57 : theA(0), theZ(0), theL(0), aEff(0.0), zEff << 48 : theA(0), theZ(0), aEff(0.0), zEff(0) 58 { 49 { 59 pnBlackTrackEnergy = 0.0; 50 pnBlackTrackEnergy = 0.0; 60 dtaBlackTrackEnergy = 0.0; 51 dtaBlackTrackEnergy = 0.0; 61 pnBlackTrackEnergyfromAnnihilation = 0.0; 52 pnBlackTrackEnergyfromAnnihilation = 0.0; 62 dtaBlackTrackEnergyfromAnnihilation = 0.0; 53 dtaBlackTrackEnergyfromAnnihilation = 0.0; 63 excitationEnergy = 0.0; 54 excitationEnergy = 0.0; 64 momentum = G4ThreeVector(0.,0.,0.); 55 momentum = G4ThreeVector(0.,0.,0.); 65 fermiMomentum = 1.52*hbarc/fermi; 56 fermiMomentum = 1.52*hbarc/fermi; 66 theTemp = 293.16*kelvin; 57 theTemp = 293.16*kelvin; 67 fIsotope = 0; << 68 } 58 } 69 59 70 G4Nucleus::G4Nucleus( const G4double A, const << 60 G4Nucleus::G4Nucleus( const G4double A, const G4double Z ) 71 { 61 { 72 SetParameters( A, Z, std::max(numberOfLambda << 62 SetParameters( A, Z ); 73 pnBlackTrackEnergy = 0.0; 63 pnBlackTrackEnergy = 0.0; 74 dtaBlackTrackEnergy = 0.0; 64 dtaBlackTrackEnergy = 0.0; 75 pnBlackTrackEnergyfromAnnihilation = 0.0; 65 pnBlackTrackEnergyfromAnnihilation = 0.0; 76 dtaBlackTrackEnergyfromAnnihilation = 0.0; 66 dtaBlackTrackEnergyfromAnnihilation = 0.0; 77 excitationEnergy = 0.0; 67 excitationEnergy = 0.0; 78 momentum = G4ThreeVector(0.,0.,0.); 68 momentum = G4ThreeVector(0.,0.,0.); 79 fermiMomentum = 1.52*hbarc/fermi; 69 fermiMomentum = 1.52*hbarc/fermi; 80 theTemp = 293.16*kelvin; 70 theTemp = 293.16*kelvin; 81 fIsotope = 0; << 82 } 71 } 83 72 84 G4Nucleus::G4Nucleus( const G4int A, const G4i << 73 G4Nucleus::G4Nucleus( const G4int A, const G4int Z ) 85 { 74 { 86 SetParameters( A, Z, std::max(numberOfLambda << 75 SetParameters( A, Z ); 87 pnBlackTrackEnergy = 0.0; 76 pnBlackTrackEnergy = 0.0; 88 dtaBlackTrackEnergy = 0.0; 77 dtaBlackTrackEnergy = 0.0; 89 pnBlackTrackEnergyfromAnnihilation = 0.0; 78 pnBlackTrackEnergyfromAnnihilation = 0.0; 90 dtaBlackTrackEnergyfromAnnihilation = 0.0; 79 dtaBlackTrackEnergyfromAnnihilation = 0.0; 91 excitationEnergy = 0.0; 80 excitationEnergy = 0.0; 92 momentum = G4ThreeVector(0.,0.,0.); 81 momentum = G4ThreeVector(0.,0.,0.); 93 fermiMomentum = 1.52*hbarc/fermi; 82 fermiMomentum = 1.52*hbarc/fermi; 94 theTemp = 293.16*kelvin; 83 theTemp = 293.16*kelvin; 95 fIsotope = 0; << 96 } 84 } 97 85 98 G4Nucleus::G4Nucleus( const G4Material *aMater 86 G4Nucleus::G4Nucleus( const G4Material *aMaterial ) 99 { 87 { 100 ChooseParameters( aMaterial ); 88 ChooseParameters( aMaterial ); 101 pnBlackTrackEnergy = 0.0; 89 pnBlackTrackEnergy = 0.0; 102 dtaBlackTrackEnergy = 0.0; 90 dtaBlackTrackEnergy = 0.0; 103 pnBlackTrackEnergyfromAnnihilation = 0.0; 91 pnBlackTrackEnergyfromAnnihilation = 0.0; 104 dtaBlackTrackEnergyfromAnnihilation = 0.0; 92 dtaBlackTrackEnergyfromAnnihilation = 0.0; 105 excitationEnergy = 0.0; 93 excitationEnergy = 0.0; 106 momentum = G4ThreeVector(0.,0.,0.); 94 momentum = G4ThreeVector(0.,0.,0.); 107 fermiMomentum = 1.52*hbarc/fermi; 95 fermiMomentum = 1.52*hbarc/fermi; 108 theTemp = aMaterial->GetTemperature(); 96 theTemp = aMaterial->GetTemperature(); 109 fIsotope = 0; << 110 } 97 } 111 98 112 G4Nucleus::~G4Nucleus() {} 99 G4Nucleus::~G4Nucleus() {} 113 100 114 << 101 G4ReactionProduct G4Nucleus:: 115 //-------------------------------------------- << 102 GetBiasedThermalNucleus(G4double aMass, G4ThreeVector aVelocity, G4double temp) const 116 // SVT (Sampling of the Velocity of the Target << 117 //-------------------------------------------- << 118 G4ReactionProduct << 119 G4Nucleus::GetBiasedThermalNucleus(G4double aM << 120 { 103 { 121 // If E_neutron <= E_threshold, Then apply t << 104 G4double velMag = aVelocity.mag(); 122 // Else consider the target nucleus being wi << 123 G4double E_threshold = G4HadronicParameters: << 124 if ( E_threshold == -1. ) { << 125 E_threshold = 400.0*8.617333262E-11*temp; << 126 } << 127 G4double E_neutron = 0.5*aVelocity.mag2()*G4 << 128 << 129 G4ReactionProduct result; 105 G4ReactionProduct result; 130 result.SetMass(aMass*G4Neutron::Neutron()->G << 106 G4double value = 0; 131 << 107 G4double random = 1; 132 if ( E_neutron <= E_threshold ) { << 108 G4double norm = 3.*std::sqrt(k_Boltzmann*temp*aMass*G4Neutron::Neutron()->GetPDGMass()); 133 << 109 norm /= G4Neutron::Neutron()->GetPDGMass(); 134 // Beta = sqrt(m/2kT) << 110 norm *= 5.; 135 G4double beta = std::sqrt(result.GetMass() << 111 norm += velMag; 136 << 112 norm /= velMag; 137 // Neutron speed vn << 113 while(value/norm<random) 138 G4double vN_norm = aVelocity.mag(); << 114 { 139 G4double vN_norm2 = vN_norm*vN_norm; << 115 result = GetThermalNucleus(aMass, temp); 140 G4double y = beta*vN_norm; << 116 G4ThreeVector targetVelocity = 1./result.GetMass()*result.GetMomentum(); 141 << 117 value = (targetVelocity+aVelocity).mag()/velMag; 142 // Normalize neutron velocity << 118 random = G4UniformRand(); 143 aVelocity = (1./vN_norm)*aVelocity; << 144 << 145 // Sample target speed << 146 G4double x2; << 147 G4double randThreshold; << 148 G4double vT_norm, vT_norm2, mu; //theta, v << 149 G4double acceptThreshold; << 150 G4double vRelativeSpeed; << 151 G4double cdf0 = 2./(2.+std::sqrt(CLHEP::pi << 152 << 153 do { << 154 // Sample the target velocity vT in the << 155 if ( G4UniformRand() < cdf0 ) { << 156 // Sample in C45 from https://laws.lan << 157 x2 = -std::log(G4UniformRand()*G4Unifo << 158 } else { << 159 // Sample in C61 from https://laws.lan << 160 G4double ampl = std::cos(CLHEP::pi/2.0 << 161 x2 = -std::log(G4UniformRand()) - std: << 162 } << 163 << 164 vT_norm = std::sqrt(x2)/beta; << 165 vT_norm2 = vT_norm*vT_norm; << 166 << 167 // Sample cosine between the incident ne << 168 mu = 2*G4UniformRand() - 1; << 169 << 170 // Define acceptance threshold << 171 vRelativeSpeed = std::sqrt(vN_norm2 + vT << 172 acceptThreshold = vRelativeSpeed/(vN_nor << 173 randThreshold = G4UniformRand(); << 174 } while ( randThreshold >= acceptThreshold << 175 << 176 DoKinematicsOfThermalNucleus(mu, vT_norm, << 177 << 178 } else { // target nucleus considered as bei << 179 << 180 result.SetMomentum(0., 0., 0.); << 181 result.SetKineticEnergy(0.); << 182 << 183 } 119 } 184 << 185 return result; 120 return result; 186 } 121 } 187 122 188 << 189 void << 190 G4Nucleus::DoKinematicsOfThermalNucleus(const << 191 G4Reac << 192 << 193 // Get target nucleus direction from the neu << 194 G4double cosTh = mu; << 195 G4ThreeVector uNorm = aVelocity; << 196 << 197 G4double sinTh = std::sqrt(1. - cosTh*cosTh) << 198 << 199 // Sample randomly the phi angle between the << 200 G4double phi = CLHEP::twopi*G4UniformRand(); << 201 G4double sinPhi = std::sin(phi); << 202 G4double cosPhi = std::cos(phi); << 203 << 204 // Find orthogonal vector to aVelocity - sol << 205 G4ThreeVector ortho(1., 1., 1.); << 206 if ( uNorm[0] ) ortho[0] = -(uNorm[1]+ << 207 else if ( uNorm[1] ) ortho[1] = -(uNorm[0]+ << 208 else if ( uNorm[2] ) ortho[2] = -(uNorm[0]+ << 209 << 210 // Normalize the vector << 211 ortho = (1/ortho.mag())*ortho; << 212 << 213 // Find vector to draw a plan perpendicular << 214 G4ThreeVector orthoComp( uNorm[1]*ortho[2] - << 215 uNorm[2]*ortho[0] - << 216 uNorm[0]*ortho[1] - << 217 << 218 // Find the direction of the target velocity << 219 G4ThreeVector directionTarget( cosTh*uNorm[0 << 220 cosTh*uNorm[1 << 221 cosTh*uNorm[2 << 222 << 223 // Normalize directionTarget << 224 directionTarget = ( 1./directionTarget.mag() << 225 << 226 // Set momentum << 227 G4double px = result.GetMass()*vT_norm*direc << 228 G4double py = result.GetMass()*vT_norm*direc << 229 G4double pz = result.GetMass()*vT_norm*direc << 230 result.SetMomentum(px, py, pz); << 231 << 232 G4double tMom = std::sqrt(px*px+py*py+pz*pz) << 233 G4double tEtot = std::sqrt( (tMom+result.Get << 234 - 2.*tMom*result.GetMass() << 235 << 236 if ( tEtot/result.GetMass() - 1. > 0.001 ) { << 237 // use relativistic energy for higher ener << 238 result.SetTotalEnergy(tEtot); << 239 } else { << 240 // use p**2/2M for lower energies (to pres << 241 result.SetKineticEnergy(tMom*tMom/(2.*resu << 242 } << 243 << 244 } << 245 << 246 << 247 G4ReactionProduct 123 G4ReactionProduct 248 G4Nucleus::GetThermalNucleus(G4double targetMa 124 G4Nucleus::GetThermalNucleus(G4double targetMass, G4double temp) const 249 { 125 { 250 G4double currentTemp = temp; << 126 G4double currentTemp = temp; 251 if (currentTemp < 0) currentTemp = theTemp; << 127 if(currentTemp < 0) currentTemp = theTemp; 252 G4ReactionProduct theTarget; << 128 G4ReactionProduct theTarget; 253 theTarget.SetMass(targetMass*G4Neutron::Neut << 129 theTarget.SetMass(targetMass*G4Neutron::Neutron()->GetPDGMass()); 254 G4double px, py, pz; << 130 G4double px, py, pz; 255 px = GetThermalPz(theTarget.GetMass(), curre << 131 px = GetThermalPz(theTarget.GetMass(), currentTemp); 256 py = GetThermalPz(theTarget.GetMass(), curre << 132 py = GetThermalPz(theTarget.GetMass(), currentTemp); 257 pz = GetThermalPz(theTarget.GetMass(), curre << 133 pz = GetThermalPz(theTarget.GetMass(), currentTemp); 258 theTarget.SetMomentum(px, py, pz); << 134 theTarget.SetMomentum(px, py, pz); 259 G4double tMom = std::sqrt(px*px+py*py+pz*pz) << 135 G4double tMom = std::sqrt(px*px+py*py+pz*pz); 260 G4double tEtot = std::sqrt((tMom+theTarget.G << 136 G4double tEtot = std::sqrt((tMom+theTarget.GetMass())* 261 (tMom+theTarget.G << 137 (tMom+theTarget.GetMass())- 262 2.*tMom*theTarge << 138 2.*tMom*theTarget.GetMass()); 263 // if(1-tEtot/theTarget.GetMass()>0.001) t << 139 if(1-tEtot/theTarget.GetMass()>0.001) 264 if (tEtot/theTarget.GetMass() - 1. > 0.001) << 140 { 265 // use relativistic energy for higher ener << 141 theTarget.SetTotalEnergy(tEtot); 266 theTarget.SetTotalEnergy(tEtot); << 142 } 267 << 143 else 268 } else { << 144 { 269 // use p**2/2M for lower energies (to pres << 145 theTarget.SetKineticEnergy(tMom*tMom/(2.*theTarget.GetMass())); 270 theTarget.SetKineticEnergy(tMom*tMom/(2.*t << 146 } 271 } << 147 return theTarget; 272 return theTarget; << 273 } 148 } 274 149 275 150 276 void 151 void 277 G4Nucleus::ChooseParameters(const G4Material* 152 G4Nucleus::ChooseParameters(const G4Material* aMaterial) 278 { 153 { 279 G4double random = G4UniformRand(); 154 G4double random = G4UniformRand(); 280 G4double sum = aMaterial->GetTotNbOfAtomsPer 155 G4double sum = aMaterial->GetTotNbOfAtomsPerVolume(); 281 const G4ElementVector* theElementVector = aM 156 const G4ElementVector* theElementVector = aMaterial->GetElementVector(); 282 G4double running(0); 157 G4double running(0); 283 // G4Element* element(0); 158 // G4Element* element(0); 284 const G4Element* element = (*theElementVecto << 159 G4Element* element = (*theElementVector)[aMaterial->GetNumberOfElements()-1]; 285 160 286 for (unsigned int i = 0; i < aMaterial->GetN 161 for (unsigned int i = 0; i < aMaterial->GetNumberOfElements(); ++i) { 287 running += aMaterial->GetVecNbOfAtomsPerVo 162 running += aMaterial->GetVecNbOfAtomsPerVolume()[i]; 288 if (running > random*sum) { 163 if (running > random*sum) { 289 element = (*theElementVector)[i]; 164 element = (*theElementVector)[i]; 290 break; 165 break; 291 } 166 } 292 } 167 } 293 168 294 if (element->GetNumberOfIsotopes() > 0) { 169 if (element->GetNumberOfIsotopes() > 0) { 295 G4double randomAbundance = G4UniformRand() 170 G4double randomAbundance = G4UniformRand(); 296 G4double sumAbundance = element->GetRelati 171 G4double sumAbundance = element->GetRelativeAbundanceVector()[0]; 297 unsigned int iso=0; 172 unsigned int iso=0; 298 while (iso < element->GetNumberOfIsotopes( << 173 while (iso < element->GetNumberOfIsotopes() && 299 sumAbundance < randomAbundance) { 174 sumAbundance < randomAbundance) { 300 ++iso; 175 ++iso; 301 sumAbundance += element->GetRelativeAbun 176 sumAbundance += element->GetRelativeAbundanceVector()[iso]; 302 } 177 } 303 theA=element->GetIsotope(iso)->GetN(); 178 theA=element->GetIsotope(iso)->GetN(); 304 theZ=element->GetIsotope(iso)->GetZ(); 179 theZ=element->GetIsotope(iso)->GetZ(); 305 theL=0; << 306 aEff=theA; 180 aEff=theA; 307 zEff=theZ; 181 zEff=theZ; 308 } else { 182 } else { 309 aEff = element->GetN(); 183 aEff = element->GetN(); 310 zEff = element->GetZ(); 184 zEff = element->GetZ(); 311 theZ = G4int(zEff + 0.5); 185 theZ = G4int(zEff + 0.5); 312 theA = G4int(aEff + 0.5); << 186 theA = G4int(aEff + 0.5); 313 theL=0; << 314 } 187 } 315 } 188 } 316 189 317 190 318 void 191 void 319 G4Nucleus::SetParameters( const G4double A, co << 192 G4Nucleus::SetParameters(const G4double A, const G4double Z) 320 { 193 { 321 theZ = G4lrint(Z); << 194 theZ = G4int(Z + 0.5); 322 theA = G4lrint(A); << 195 theA = G4int(A + 0.5); 323 theL = std::max(numberOfLambdas, 0); << 324 if (theA<1 || theZ<0 || theZ>theA) { 196 if (theA<1 || theZ<0 || theZ>theA) { 325 throw G4HadronicException(__FILE__, __LINE 197 throw G4HadronicException(__FILE__, __LINE__, 326 "G4Nucleus::SetParameters called w 198 "G4Nucleus::SetParameters called with non-physical parameters"); 327 } 199 } 328 aEff = A; // atomic weight 200 aEff = A; // atomic weight 329 zEff = Z; // atomic number 201 zEff = Z; // atomic number 330 fIsotope = 0; << 331 } 202 } 332 203 333 << 204 void 334 void << 205 G4Nucleus::SetParameters( const G4int A, const G4int Z ) 335 G4Nucleus::SetParameters( const G4int A, const << 206 { 336 { << 207 theZ = Z; 337 theZ = Z; << 208 theA = A; 338 theA = A; << 209 if( theA<1 || theZ<0 || theZ>theA ) 339 theL = std::max(numberOfLambdas, 0); << 340 if( theA<1 || theZ<0 || theZ>theA ) << 341 { 210 { 342 throw G4HadronicException(__FILE__, __LI 211 throw G4HadronicException(__FILE__, __LINE__, 343 "G4Nucleus::SetParameters called with << 212 "G4Nucleus::SetParameters called with non-physical parameters"); 344 } 213 } 345 aEff = A; // atomic weight << 214 aEff = A; // atomic weight 346 zEff = Z; // atomic number << 215 zEff = Z; // atomic number 347 fIsotope = 0; << 348 } << 349 << 350 << 351 G4DynamicParticle * << 352 G4Nucleus::ReturnTargetParticle() const << 353 { << 354 // choose a proton or a neutron (or a lamba << 355 G4DynamicParticle *targetParticle = new G4Dy << 356 const G4double rnd = G4UniformRand(); << 357 if ( rnd < zEff/aEff ) { << 358 targetParticle->SetDefinition( G4Proton::P << 359 } else if ( rnd < (zEff + theL*1.0)/aEff ) { << 360 targetParticle->SetDefinition( G4Lambda::L << 361 } else { << 362 targetParticle->SetDefinition( G4Neutron:: << 363 } 216 } 364 return targetParticle; << 365 } << 366 217 >> 218 G4DynamicParticle * >> 219 G4Nucleus::ReturnTargetParticle() const >> 220 { >> 221 // choose a proton or a neutron as the target particle >> 222 >> 223 G4DynamicParticle *targetParticle = new G4DynamicParticle; >> 224 if( G4UniformRand() < zEff/aEff ) >> 225 targetParticle->SetDefinition( G4Proton::Proton() ); >> 226 else >> 227 targetParticle->SetDefinition( G4Neutron::Neutron() ); >> 228 return targetParticle; >> 229 } 367 230 368 G4double << 231 G4double 369 G4Nucleus::AtomicMass( const G4double A, const << 232 G4Nucleus::AtomicMass( const G4double A, const G4double Z ) const 370 { << 233 { 371 // Now returns (atomic mass - electron masse << 234 // Now returns (atomic mass - electron masses) 372 if ( numberOfLambdas > 0 ) { << 373 return G4HyperNucleiProperties::GetNuclear << 374 } else { << 375 return G4NucleiProperties::GetNuclearMass( 235 return G4NucleiProperties::GetNuclearMass(A, Z); 376 } 236 } 377 } << 378 << 379 237 380 G4double << 238 G4double 381 G4Nucleus::AtomicMass( const G4int A, const G4 << 239 G4Nucleus::AtomicMass( const G4int A, const G4int Z ) const 382 { << 240 { 383 // Now returns (atomic mass - electron masse << 241 // Now returns (atomic mass - electron masses) 384 if ( numberOfLambdas > 0 ) { << 385 return G4HyperNucleiProperties::GetNuclear << 386 } else { << 387 return G4NucleiProperties::GetNuclearMass( 242 return G4NucleiProperties::GetNuclearMass(A, Z); 388 } 243 } 389 } << 390 244 >> 245 G4double >> 246 G4Nucleus::GetThermalPz( const G4double mass, const G4double temp ) const >> 247 { >> 248 G4double result = G4RandGauss::shoot(); >> 249 result *= std::sqrt(k_Boltzmann*temp*mass); // Das ist impuls (Pz), >> 250 // nichtrelativistische rechnung >> 251 // Maxwell verteilung angenommen >> 252 return result; >> 253 } 391 254 392 G4double << 255 G4double 393 G4Nucleus::GetThermalPz( const G4double mass, << 256 G4Nucleus::EvaporationEffects( G4double kineticEnergy ) 394 { << 257 { 395 G4double result = G4RandGauss::shoot(); << 258 // derived from original FORTRAN code EXNU by H. Fesefeldt (10-Dec-1986) 396 result *= std::sqrt(k_Boltzmann*temp*mass); << 259 // 397 // ni << 260 // Nuclear evaporation as function of atomic number 398 // Ma << 261 // and kinetic energy (MeV) of primary particle 399 return result; << 262 // 400 } << 263 // returns kinetic energy (MeV) >> 264 // >> 265 if( aEff < 1.5 ) >> 266 { >> 267 pnBlackTrackEnergy = dtaBlackTrackEnergy = 0.0; >> 268 return 0.0; >> 269 } >> 270 G4double ek = kineticEnergy/GeV; >> 271 G4float ekin = std::min( 4.0, std::max( 0.1, ek ) ); >> 272 const G4float atno = std::min( 120., aEff ); >> 273 const G4float gfa = 2.0*((aEff-1.0)/70.)*std::exp(-(aEff-1.0)/70.); >> 274 // >> 275 // 0.35 value at 1 GeV >> 276 // 0.05 value at 0.1 GeV >> 277 // >> 278 G4float cfa = std::max( 0.15, 0.35 + ((0.35-0.05)/2.3)*std::log(ekin) ); >> 279 G4float exnu = 7.716 * cfa * std::exp(-cfa) >> 280 * ((atno-1.0)/120.)*std::exp(-(atno-1.0)/120.); >> 281 G4float fpdiv = std::max( 0.5, 1.0-0.25*ekin*ekin ); >> 282 // >> 283 // pnBlackTrackEnergy is the kinetic energy (in GeV) available for >> 284 // proton/neutron black track particles >> 285 // dtaBlackTrackEnergy is the kinetic energy (in GeV) available for >> 286 // deuteron/triton/alpha black track particles >> 287 // >> 288 pnBlackTrackEnergy = exnu*fpdiv; >> 289 dtaBlackTrackEnergy = exnu*(1.0-fpdiv); >> 290 >> 291 if( G4int(zEff+0.1) != 82 ) >> 292 { >> 293 G4double ran1 = -6.0; >> 294 G4double ran2 = -6.0; >> 295 for( G4int i=0; i<12; ++i ) >> 296 { >> 297 ran1 += G4UniformRand(); >> 298 ran2 += G4UniformRand(); >> 299 } >> 300 pnBlackTrackEnergy *= 1.0 + ran1*gfa; >> 301 dtaBlackTrackEnergy *= 1.0 + ran2*gfa; >> 302 } >> 303 pnBlackTrackEnergy = std::max( 0.0, pnBlackTrackEnergy ); >> 304 dtaBlackTrackEnergy = std::max( 0.0, dtaBlackTrackEnergy ); >> 305 while( pnBlackTrackEnergy+dtaBlackTrackEnergy >= ek ) >> 306 { >> 307 pnBlackTrackEnergy *= 1.0 - 0.5*G4UniformRand(); >> 308 dtaBlackTrackEnergy *= 1.0 - 0.5*G4UniformRand(); >> 309 } >> 310 // G4cout << "EvaporationEffects "<<kineticEnergy<<" " >> 311 // <<pnBlackTrackEnergy+dtaBlackTrackEnergy<<endl; >> 312 return (pnBlackTrackEnergy+dtaBlackTrackEnergy)*GeV; >> 313 } 401 314 >> 315 G4double G4Nucleus::AnnihilationEvaporationEffects(G4double kineticEnergy, G4double ekOrg) >> 316 { >> 317 // Nuclear evaporation as a function of atomic number and kinetic >> 318 // energy (MeV) of primary particle. Modified for annihilation effects. >> 319 // >> 320 if( aEff < 1.5 || ekOrg < 0.) >> 321 { >> 322 pnBlackTrackEnergyfromAnnihilation = 0.0; >> 323 dtaBlackTrackEnergyfromAnnihilation = 0.0; >> 324 return 0.0; >> 325 } >> 326 G4double ek = kineticEnergy/GeV; >> 327 G4float ekin = std::min( 4.0, std::max( 0.1, ek ) ); >> 328 const G4float atno = std::min( 120., aEff ); >> 329 const G4float gfa = 2.0*((aEff-1.0)/70.)*std::exp(-(aEff-1.0)/70.); >> 330 >> 331 G4float cfa = std::max( 0.15, 0.35 + ((0.35-0.05)/2.3)*std::log(ekin) ); >> 332 G4float exnu = 7.716 * cfa * std::exp(-cfa) >> 333 * ((atno-1.0)/120.)*std::exp(-(atno-1.0)/120.); >> 334 G4float fpdiv = std::max( 0.5, 1.0-0.25*ekin*ekin ); 402 335 403 G4double << 336 pnBlackTrackEnergyfromAnnihilation = exnu*fpdiv; 404 G4Nucleus::EvaporationEffects( G4double kineti << 337 dtaBlackTrackEnergyfromAnnihilation = exnu*(1.0-fpdiv); 405 { << 338 406 // derived from original FORTRAN code EXNU b << 407 // << 408 // Nuclear evaporation as function of atomic << 409 // and kinetic energy (MeV) of primary parti << 410 // << 411 // returns kinetic energy (MeV) << 412 // << 413 if( aEff < 1.5 ) << 414 { << 415 pnBlackTrackEnergy = dtaBlackTrackEnergy = << 416 return 0.0; << 417 } << 418 G4double ek = kineticEnergy/GeV; << 419 G4float ekin = std::min( 4.0, std::max( 0.1, << 420 const G4float atno = std::min( 120., aEff ); << 421 const G4float gfa = 2.0*((aEff-1.0)/70.)*G4E << 422 // << 423 // 0.35 value at 1 GeV << 424 // 0.05 value at 0.1 GeV << 425 // << 426 G4float cfa = std::max( 0.15, 0.35 + ((0.35- << 427 G4float exnu = 7.716 * cfa * G4Exp(-cfa) << 428 * ((atno-1.0)/120.)*G4Exp(-(atno-1.0)/120. << 429 G4float fpdiv = std::max( 0.5, 1.0-0.25*ekin << 430 // << 431 // pnBlackTrackEnergy is the kinetic energy << 432 // proton/neutron black << 433 // dtaBlackTrackEnergy is the kinetic energy << 434 // deuteron/triton/alpha << 435 // << 436 pnBlackTrackEnergy = exnu*fpdiv; << 437 dtaBlackTrackEnergy = exnu*(1.0-fpdiv); << 438 << 439 if( G4int(zEff+0.1) != 82 ) << 440 { << 441 G4double ran1 = -6.0; 339 G4double ran1 = -6.0; 442 G4double ran2 = -6.0; 340 G4double ran2 = -6.0; 443 for( G4int i=0; i<12; ++i ) << 341 for( G4int i=0; i<12; ++i ) { 444 { << 445 ran1 += G4UniformRand(); 342 ran1 += G4UniformRand(); 446 ran2 += G4UniformRand(); 343 ran2 += G4UniformRand(); 447 } 344 } 448 pnBlackTrackEnergy *= 1.0 + ran1*gfa; << 345 pnBlackTrackEnergyfromAnnihilation *= 1.0 + ran1*gfa; 449 dtaBlackTrackEnergy *= 1.0 + ran2*gfa; << 346 dtaBlackTrackEnergyfromAnnihilation *= 1.0 + ran2*gfa; 450 } << 451 pnBlackTrackEnergy = std::max( 0.0, pnBlackT << 452 dtaBlackTrackEnergy = std::max( 0.0, dtaBlac << 453 while( pnBlackTrackEnergy+dtaBlackTrackEnerg << 454 { << 455 pnBlackTrackEnergy *= 1.0 - 0.5*G4UniformR << 456 dtaBlackTrackEnergy *= 1.0 - 0.5*G4Uniform << 457 } << 458 //G4cout << "EvaporationEffects "<<kineticEn << 459 // <<pnBlackTrackEnergy+dtaBlackTrackE << 460 return (pnBlackTrackEnergy+dtaBlackTrackEner << 461 } << 462 << 463 << 464 G4double << 465 G4Nucleus::AnnihilationEvaporationEffects(G4do << 466 { << 467 // Nuclear evaporation as a function of atom << 468 // energy (MeV) of primary particle. Modifi << 469 // << 470 if( aEff < 1.5 || ekOrg < 0.) << 471 { << 472 pnBlackTrackEnergyfromAnnihilation = 0.0; << 473 dtaBlackTrackEnergyfromAnnihilation = 0.0; << 474 return 0.0; << 475 } << 476 G4double ek = kineticEnergy/GeV; << 477 G4float ekin = std::min( 4.0, std::max( 0.1, << 478 const G4float atno = std::min( 120., aEff ); << 479 const G4float gfa = 2.0*((aEff-1.0)/70.)*G4E << 480 << 481 G4float cfa = std::max( 0.15, 0.35 + ((0.35- << 482 G4float exnu = 7.716 * cfa * G4Exp(-cfa) << 483 * ((atno-1.0)/120.)*G4Exp(-(atno-1.0)/120. << 484 G4float fpdiv = std::max( 0.5, 1.0-0.25*ekin << 485 << 486 pnBlackTrackEnergyfromAnnihilation = exnu*fp << 487 dtaBlackTrackEnergyfromAnnihilation = exnu*( << 488 << 489 G4double ran1 = -6.0; << 490 G4double ran2 = -6.0; << 491 for( G4int i=0; i<12; ++i ) { << 492 ran1 += G4UniformRand(); << 493 ran2 += G4UniformRand(); << 494 } << 495 pnBlackTrackEnergyfromAnnihilation *= 1.0 + << 496 dtaBlackTrackEnergyfromAnnihilation *= 1.0 + << 497 << 498 pnBlackTrackEnergyfromAnnihilation = std::ma << 499 dtaBlackTrackEnergyfromAnnihilation = std::m << 500 G4double blackSum = pnBlackTrackEnergyfromAn << 501 if (blackSum >= ekOrg/GeV) { << 502 pnBlackTrackEnergyfromAnnihilation *= ekOr << 503 dtaBlackTrackEnergyfromAnnihilation *= ekO << 504 } << 505 347 506 return (pnBlackTrackEnergyfromAnnihilation+d << 348 pnBlackTrackEnergyfromAnnihilation = std::max( 0.0, pnBlackTrackEnergyfromAnnihilation); 507 } << 349 dtaBlackTrackEnergyfromAnnihilation = std::max( 0.0, dtaBlackTrackEnergyfromAnnihilation); >> 350 G4double blackSum = pnBlackTrackEnergyfromAnnihilation+dtaBlackTrackEnergyfromAnnihilation; >> 351 if (blackSum >= ekOrg/GeV) { >> 352 pnBlackTrackEnergyfromAnnihilation *= ekOrg/GeV/blackSum; >> 353 dtaBlackTrackEnergyfromAnnihilation *= ekOrg/GeV/blackSum; >> 354 } 508 355 >> 356 return (pnBlackTrackEnergyfromAnnihilation+dtaBlackTrackEnergyfromAnnihilation)*GeV; >> 357 } 509 358 510 G4double << 359 G4double 511 G4Nucleus::Cinema( G4double kineticEnergy ) << 360 G4Nucleus::Cinema( G4double kineticEnergy ) 512 { << 361 { 513 // derived from original FORTRAN code CINEMA << 362 // derived from original FORTRAN code CINEMA by H. Fesefeldt (14-Oct-1987) 514 // << 363 // 515 // input: kineticEnergy (MeV) << 364 // input: kineticEnergy (MeV) 516 // returns modified kinetic energy (MeV) << 365 // returns modified kinetic energy (MeV) 517 // << 366 // 518 static const G4double expxu = 82.; << 367 static const G4double expxu = 82.; // upper bound for arg. of exp 519 static const G4double expxl = -expxu; << 368 static const G4double expxl = -expxu; // lower bound for arg. of exp 520 << 369 521 G4double ek = kineticEnergy/GeV; << 370 G4double ek = kineticEnergy/GeV; 522 G4double ekLog = G4Log( ek ); << 371 G4double ekLog = std::log( ek ); 523 G4double aLog = G4Log( aEff ); << 372 G4double aLog = std::log( aEff ); 524 G4double em = std::min( 1.0, 0.2390 + 0.0408 << 373 G4double em = std::min( 1.0, 0.2390 + 0.0408*aLog*aLog ); 525 G4double temp1 = -ek * std::min( 0.15, 0.001 << 374 G4double temp1 = -ek * std::min( 0.15, 0.0019*aLog*aLog*aLog ); 526 G4double temp2 = G4Exp( std::max( expxl, std << 375 G4double temp2 = std::exp( std::max( expxl, std::min( expxu, -(ekLog-em)*(ekLog-em)*2.0 ) ) ); 527 G4double result = 0.0; << 376 G4double result = 0.0; 528 if( std::abs( temp1 ) < 1.0 ) << 377 if( std::abs( temp1 ) < 1.0 ) 529 { << 378 { 530 if( temp2 > 1.0e-10 )result = temp1*temp2; << 379 if( temp2 > 1.0e-10 )result = temp1*temp2; 531 } << 380 } 532 else result = temp1*temp2; << 381 else result = temp1*temp2; 533 if( result < -ek )result = -ek; << 382 if( result < -ek )result = -ek; 534 return result*GeV; << 383 return result*GeV; 535 } << 384 } 536 385 >> 386 // >> 387 // methods for class G4Nucleus ... by Christian Volcker >> 388 // 537 389 538 G4ThreeVector G4Nucleus::GetFermiMomentum() << 390 G4ThreeVector G4Nucleus::GetFermiMomentum() 539 { << 391 { 540 // chv: .. we assume zero temperature! << 392 // chv: .. we assume zero temperature! 541 << 393 542 // momentum is equally distributed in each p << 394 // momentum is equally distributed in each phasespace volume dpx, dpy, dpz. 543 G4double ranflat1= << 395 G4double ranflat1= 544 G4RandFlat::shoot((G4double)0.,(G4double)f << 396 CLHEP::RandFlat::shoot((G4double)0.,(G4double)fermiMomentum); 545 G4double ranflat2= << 397 G4double ranflat2= 546 G4RandFlat::shoot((G4double)0.,(G4double)f << 398 CLHEP::RandFlat::shoot((G4double)0.,(G4double)fermiMomentum); 547 G4double ranflat3= << 399 G4double ranflat3= 548 G4RandFlat::shoot((G4double)0.,(G4double)f << 400 CLHEP::RandFlat::shoot((G4double)0.,(G4double)fermiMomentum); 549 G4double ranmax = (ranflat1>ranflat2? ranfla << 401 G4double ranmax = (ranflat1>ranflat2? ranflat1: ranflat2); 550 ranmax = (ranmax>ranflat3? ranmax : ranflat3 << 402 ranmax = (ranmax>ranflat3? ranmax : ranflat3); 551 << 403 552 // Isotropic momentum distribution << 404 // Isotropic momentum distribution 553 G4double costheta = 2.*G4UniformRand() - 1.0 << 405 G4double costheta = 2.*G4UniformRand() - 1.0; 554 G4double sintheta = std::sqrt(1.0 - costheta << 406 G4double sintheta = std::sqrt(1.0 - costheta*costheta); 555 G4double phi = 2.0*pi*G4UniformRand(); << 407 G4double phi = 2.0*pi*G4UniformRand(); 556 << 408 557 G4double pz=costheta*ranmax; << 409 G4double pz=costheta*ranmax; 558 G4double px=sintheta*std::cos(phi)*ranmax; << 410 G4double px=sintheta*std::cos(phi)*ranmax; 559 G4double py=sintheta*std::sin(phi)*ranmax; << 411 G4double py=sintheta*std::sin(phi)*ranmax; 560 G4ThreeVector p(px,py,pz); << 412 G4ThreeVector p(px,py,pz); 561 return p; << 413 return p; 562 } << 414 } 563 415 564 << 416 G4ReactionProductVector* G4Nucleus::Fragmentate() 565 G4ReactionProductVector* G4Nucleus::Fragmentat << 417 { 566 { << 418 // needs implementation! 567 // needs implementation! << 419 return NULL; 568 return nullptr; << 420 } 569 } << 570 421 571 << 422 void G4Nucleus::AddMomentum(const G4ThreeVector aMomentum) 572 void G4Nucleus::AddMomentum(const G4ThreeVecto << 423 { 573 { << 424 momentum+=(aMomentum); 574 momentum+=(aMomentum); << 425 } 575 } << 576 << 577 426 578 void G4Nucleus::AddExcitationEnergy( G4double << 427 void G4Nucleus::AddExcitationEnergy( G4double anEnergy ) 579 { << 428 { 580 excitationEnergy+=anEnergy; << 429 excitationEnergy+=anEnergy; 581 } << 430 } 582 431 583 /* end of file */ 432 /* end of file */ 584 433 585 434