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 // G4 Model: Low Energy Fission 28 // G4 Model: Low Energy Fission 29 // F.W. Jones, TRIUMF, 03-DEC-96 29 // F.W. Jones, TRIUMF, 03-DEC-96 30 // 30 // 31 // This is a prototype of a low-energy fission 31 // This is a prototype of a low-energy fission process. 32 // Currently it is based on the GHEISHA routin 32 // Currently it is based on the GHEISHA routine FISSIO, 33 // and conforms fairly closely to the original 33 // and conforms fairly closely to the original Fortran. 34 // Note: energy is in MeV and momentum is in M 34 // Note: energy is in MeV and momentum is in MeV/c. 35 // 35 // 36 // use -scheme for elastic scattering: HPW, 20 36 // use -scheme for elastic scattering: HPW, 20th June 1997 37 // the code comes mostly from the old Low-ener 37 // the code comes mostly from the old Low-energy Fission class 38 // 38 // 39 // 25-JUN-98 FWJ: replaced missing Initialize 39 // 25-JUN-98 FWJ: replaced missing Initialize for ParticleChange. 40 40 41 #include <iostream> 41 #include <iostream> 42 42 43 #include "G4LFission.hh" 43 #include "G4LFission.hh" 44 #include "globals.hh" 44 #include "globals.hh" 45 #include "G4Exp.hh" 45 #include "G4Exp.hh" 46 #include "G4Log.hh" 46 #include "G4Log.hh" 47 #include "G4Pow.hh" 47 #include "G4Pow.hh" 48 #include "G4PhysicalConstants.hh" 48 #include "G4PhysicalConstants.hh" 49 #include "G4SystemOfUnits.hh" 49 #include "G4SystemOfUnits.hh" 50 #include "Randomize.hh" 50 #include "Randomize.hh" 51 #include "G4PhysicsModelCatalog.hh" 51 #include "G4PhysicsModelCatalog.hh" 52 52 53 G4LFission::G4LFission(const G4String& name) 53 G4LFission::G4LFission(const G4String& name) 54 : G4HadronicInteraction(name), secID(-1) 54 : G4HadronicInteraction(name), secID(-1) 55 { 55 { 56 init(); 56 init(); 57 SetMinEnergy(0.0*GeV); 57 SetMinEnergy(0.0*GeV); 58 SetMaxEnergy(DBL_MAX); 58 SetMaxEnergy(DBL_MAX); 59 G4PhysicsModelCatalog::GetModelID( "model_" 59 G4PhysicsModelCatalog::GetModelID( "model_" + GetModelName() ); 60 } 60 } 61 61 62 62 63 G4LFission::~G4LFission() 63 G4LFission::~G4LFission() 64 { 64 { 65 theParticleChange.Clear(); 65 theParticleChange.Clear(); 66 } 66 } 67 67 68 68 69 void G4LFission::ModelDescription(std::ostream 69 void G4LFission::ModelDescription(std::ostream& outFile) const 70 { 70 { 71 outFile << "G4LFission is one of the Low Ene 71 outFile << "G4LFission is one of the Low Energy Parameterized\n" 72 << "(LEP) models used to implement n 72 << "(LEP) models used to implement neutron-induced fission of\n" 73 << "nuclei. It is a re-engineered v 73 << "nuclei. It is a re-engineered version of the GHEISHA code\n" 74 << "of H. Fesefeldt which emits neut 74 << "of H. Fesefeldt which emits neutrons and gammas but no\n" 75 << "nuclear fragments. The model is 75 << "nuclear fragments. The model is applicable to all incident\n" 76 << "neutron energies.\n"; 76 << "neutron energies.\n"; 77 } 77 } 78 78 79 void G4LFission::init() 79 void G4LFission::init() 80 { 80 { 81 G4int i; 81 G4int i; 82 G4double xx = 1. - 0.5; 82 G4double xx = 1. - 0.5; 83 G4double xxx = std::sqrt(2.29*xx); 83 G4double xxx = std::sqrt(2.29*xx); 84 spneut[0] = G4Exp(-xx/0.965)*(G4Exp(xxx) - 84 spneut[0] = G4Exp(-xx/0.965)*(G4Exp(xxx) - G4Exp(-xxx))/2.; 85 for (i = 2; i <= 10; i++) { 85 for (i = 2; i <= 10; i++) { 86 xx = i*1. - 0.5; 86 xx = i*1. - 0.5; 87 xxx = std::sqrt(2.29*xx); 87 xxx = std::sqrt(2.29*xx); 88 spneut[i-1] = spneut[i-2] + G4Exp(-xx/0. 88 spneut[i-1] = spneut[i-2] + G4Exp(-xx/0.965)*(G4Exp(xxx) - G4Exp(-xxx))/2.; 89 } 89 } 90 for (i = 1; i <= 10; i++) { 90 for (i = 1; i <= 10; i++) { 91 spneut[i-1] = spneut[i-1]/spneut[9]; 91 spneut[i-1] = spneut[i-1]/spneut[9]; 92 if (verboseLevel > 1) G4cout << "G4LFiss 92 if (verboseLevel > 1) G4cout << "G4LFission::init: i=" << i << 93 " spneut=" << spneut[i-1] << G4endl; 93 " spneut=" << spneut[i-1] << G4endl; 94 } 94 } 95 } 95 } 96 96 97 97 98 G4HadFinalState* G4LFission::ApplyYourself(con 98 G4HadFinalState* G4LFission::ApplyYourself(const G4HadProjectile& aTrack, 99 G4N 99 G4Nucleus& targetNucleus) 100 { 100 { 101 theParticleChange.Clear(); 101 theParticleChange.Clear(); 102 const G4HadProjectile* aParticle = &aTrack; 102 const G4HadProjectile* aParticle = &aTrack; 103 103 104 G4double N = targetNucleus.GetA_asInt(); 104 G4double N = targetNucleus.GetA_asInt(); 105 G4double Z = targetNucleus.GetZ_asInt(); 105 G4double Z = targetNucleus.GetZ_asInt(); 106 theParticleChange.SetStatusChange(stopAndKil 106 theParticleChange.SetStatusChange(stopAndKill); 107 107 108 G4double P = aParticle->GetTotalMomentum()/M 108 G4double P = aParticle->GetTotalMomentum()/MeV; 109 G4double Px = aParticle->Get4Momentum().vect 109 G4double Px = aParticle->Get4Momentum().vect().x(); 110 G4double Py = aParticle->Get4Momentum().vect 110 G4double Py = aParticle->Get4Momentum().vect().y(); 111 G4double Pz = aParticle->Get4Momentum().vect 111 G4double Pz = aParticle->Get4Momentum().vect().z(); 112 G4double E = aParticle->GetTotalEnergy()/MeV 112 G4double E = aParticle->GetTotalEnergy()/MeV; 113 G4double E0 = aParticle->GetDefinition()->Ge 113 G4double E0 = aParticle->GetDefinition()->GetPDGMass()/MeV; 114 G4double Q = aParticle->GetDefinition()->Get 114 G4double Q = aParticle->GetDefinition()->GetPDGCharge(); 115 if (verboseLevel > 1) { 115 if (verboseLevel > 1) { 116 G4cout << "G4LFission:ApplyYourself: inc 116 G4cout << "G4LFission:ApplyYourself: incident particle:" << G4endl; 117 G4cout << "P " << P << " MeV/c" << 117 G4cout << "P " << P << " MeV/c" << G4endl; 118 G4cout << "Px " << Px << " MeV/c" << 118 G4cout << "Px " << Px << " MeV/c" << G4endl; 119 G4cout << "Py " << Py << " MeV/c" << 119 G4cout << "Py " << Py << " MeV/c" << G4endl; 120 G4cout << "Pz " << Pz << " MeV/c" << 120 G4cout << "Pz " << Pz << " MeV/c" << G4endl; 121 G4cout << "E " << E << " MeV" << G4 121 G4cout << "E " << E << " MeV" << G4endl; 122 G4cout << "mass " << E0 << " MeV" << G 122 G4cout << "mass " << E0 << " MeV" << G4endl; 123 G4cout << "charge " << Q << G4endl; 123 G4cout << "charge " << Q << G4endl; 124 } 124 } 125 // GHEISHA ADD operation to get total energy 125 // GHEISHA ADD operation to get total energy, mass, charge: 126 if (verboseLevel > 1) { 126 if (verboseLevel > 1) { 127 G4cout << "G4LFission:ApplyYourself: mat 127 G4cout << "G4LFission:ApplyYourself: material:" << G4endl; 128 G4cout << "A " << N << G4endl; 128 G4cout << "A " << N << G4endl; 129 G4cout << "Z " << Z << G4endl; 129 G4cout << "Z " << Z << G4endl; 130 G4cout << "atomic mass " << 130 G4cout << "atomic mass " << 131 Atomas(N, Z) << "MeV" << G4endl; 131 Atomas(N, Z) << "MeV" << G4endl; 132 } 132 } 133 E = E + Atomas(N, Z); 133 E = E + Atomas(N, Z); 134 G4double E02 = E*E - P*P; 134 G4double E02 = E*E - P*P; 135 E0 = std::sqrt(std::abs(E02)); 135 E0 = std::sqrt(std::abs(E02)); 136 if (E02 < 0) E0 = -E0; 136 if (E02 < 0) E0 = -E0; 137 Q = Q + Z; 137 Q = Q + Z; 138 if (verboseLevel > 1) { 138 if (verboseLevel > 1) { 139 G4cout << "G4LFission:ApplyYourself: tot 139 G4cout << "G4LFission:ApplyYourself: total:" << G4endl; 140 G4cout << "E " << E << " MeV" << G4 140 G4cout << "E " << E << " MeV" << G4endl; 141 G4cout << "mass " << E0 << " MeV" << G 141 G4cout << "mass " << E0 << " MeV" << G4endl; 142 G4cout << "charge " << Q << G4endl; 142 G4cout << "charge " << Q << G4endl; 143 } 143 } 144 Px = -Px; 144 Px = -Px; 145 Py = -Py; 145 Py = -Py; 146 Pz = -Pz; 146 Pz = -Pz; 147 147 148 G4double e1 = aParticle->GetKineticEnergy()/ 148 G4double e1 = aParticle->GetKineticEnergy()/MeV; 149 if (e1 < 1.) e1 = 1.; 149 if (e1 < 1.) e1 = 1.; 150 150 151 // Average number of neutrons 151 // Average number of neutrons 152 G4double avern = 2.569 + 0.559*G4Log(e1); 152 G4double avern = 2.569 + 0.559*G4Log(e1); 153 G4bool photofission = 0; // For now 153 G4bool photofission = 0; // For now 154 // Take the following value if photofission is 154 // Take the following value if photofission is not included 155 if (!photofission) avern = 2.569 + 0.900*G4 155 if (!photofission) avern = 2.569 + 0.900*G4Log(e1); 156 156 157 // Average number of gammas 157 // Average number of gammas 158 G4double averg = 9.500 + 0.600*G4Log(e1); 158 G4double averg = 9.500 + 0.600*G4Log(e1); 159 159 160 G4double ran = G4RandGauss::shoot(); 160 G4double ran = G4RandGauss::shoot(); 161 // Number of neutrons 161 // Number of neutrons 162 G4int nn = static_cast<G4int>(avern + ran*1 162 G4int nn = static_cast<G4int>(avern + ran*1.23 + 0.5); 163 ran = G4RandGauss::shoot(); 163 ran = G4RandGauss::shoot(); 164 // Number of gammas 164 // Number of gammas 165 G4int ng = static_cast<G4int>(averg + ran*3 165 G4int ng = static_cast<G4int>(averg + ran*3. + 0.5); 166 if (nn < 1) nn = 1; 166 if (nn < 1) nn = 1; 167 if (ng < 1) ng = 1; 167 if (ng < 1) ng = 1; 168 G4double exn = 0.; 168 G4double exn = 0.; 169 G4double exg = 0.; 169 G4double exg = 0.; 170 170 171 // Make secondary neutrons and distribute kine 171 // Make secondary neutrons and distribute kinetic energy 172 G4DynamicParticle* aNeutron; 172 G4DynamicParticle* aNeutron; 173 G4int i; 173 G4int i; 174 for (i = 1; i <= nn; i++) { 174 for (i = 1; i <= nn; i++) { 175 ran = G4UniformRand(); 175 ran = G4UniformRand(); 176 G4int j; 176 G4int j; 177 for (j = 1; j <= 10; j++) { 177 for (j = 1; j <= 10; j++) { 178 if (ran < spneut[j-1]) goto label12; 178 if (ran < spneut[j-1]) goto label12; 179 } 179 } 180 j = 10; 180 j = 10; 181 label12: 181 label12: 182 ran = G4UniformRand(); 182 ran = G4UniformRand(); 183 G4double ekin = (j - 1)*1. + ran; 183 G4double ekin = (j - 1)*1. + ran; 184 exn = exn + ekin; 184 exn = exn + ekin; 185 aNeutron = new G4DynamicParticle(G4Neutr 185 aNeutron = new G4DynamicParticle(G4Neutron::NeutronDefinition(), 186 G4Parti 186 G4ParticleMomentum(1.,0.,0.), 187 ekin*Me 187 ekin*MeV); 188 theParticleChange.AddSecondary(aNeutron, 188 theParticleChange.AddSecondary(aNeutron, secID); 189 } 189 } 190 190 191 // Make secondary gammas and distribute kineti 191 // Make secondary gammas and distribute kinetic energy 192 G4DynamicParticle* aGamma; 192 G4DynamicParticle* aGamma; 193 for (i = 1; i <= ng; i++) { 193 for (i = 1; i <= ng; i++) { 194 ran = G4UniformRand(); 194 ran = G4UniformRand(); 195 G4double ekin = -0.87*G4Log(ran); 195 G4double ekin = -0.87*G4Log(ran); 196 exg = exg + ekin; 196 exg = exg + ekin; 197 aGamma = new G4DynamicParticle(G4Gamma:: 197 aGamma = new G4DynamicParticle(G4Gamma::GammaDefinition(), 198 G4Particl 198 G4ParticleMomentum(1.,0.,0.), 199 ekin*MeV) 199 ekin*MeV); 200 theParticleChange.AddSecondary(aGamma, s 200 theParticleChange.AddSecondary(aGamma, secID); 201 } 201 } 202 202 203 // Distribute momentum vectors and do Lorentz 203 // Distribute momentum vectors and do Lorentz transformation 204 204 205 G4HadSecondary* theSecondary; 205 G4HadSecondary* theSecondary; 206 206 207 for (i = 1; i <= nn + ng; i++) { 207 for (i = 1; i <= nn + ng; i++) { 208 G4double ran1 = G4UniformRand(); 208 G4double ran1 = G4UniformRand(); 209 G4double ran2 = G4UniformRand(); 209 G4double ran2 = G4UniformRand(); 210 G4double cost = -1. + 2.*ran1; 210 G4double cost = -1. + 2.*ran1; 211 G4double sint = std::sqrt(std::abs(1. - 211 G4double sint = std::sqrt(std::abs(1. - cost*cost)); 212 G4double phi = ran2*twopi; 212 G4double phi = ran2*twopi; 213 // G4cout << ran1 << " " << ran2 << 213 // G4cout << ran1 << " " << ran2 << G4endl; 214 // G4cout << cost << " " << sint << 214 // G4cout << cost << " " << sint << " " << phi << G4endl; 215 theSecondary = theParticleChange.GetSeco 215 theSecondary = theParticleChange.GetSecondary(i - 1); 216 G4double pp = theSecondary->GetParticle( 216 G4double pp = theSecondary->GetParticle()->GetTotalMomentum()/MeV; 217 G4double px = pp*sint*std::sin(phi); 217 G4double px = pp*sint*std::sin(phi); 218 G4double py = pp*sint*std::cos(phi); 218 G4double py = pp*sint*std::cos(phi); 219 G4double pz = pp*cost; 219 G4double pz = pp*cost; 220 // G4cout << pp << G4endl; 220 // G4cout << pp << G4endl; 221 // G4cout << px << " " << py << " " 221 // G4cout << px << " " << py << " " << pz << G4endl; 222 G4double e = theSecondary->GetParticle() 222 G4double e = theSecondary->GetParticle()->GetTotalEnergy()/MeV; 223 G4double e0 = theSecondary->GetParticle( 223 G4double e0 = theSecondary->GetParticle()->GetDefinition()->GetPDGMass()/MeV; 224 224 225 G4double a = px*Px + py*Py + pz*Pz; 225 G4double a = px*Px + py*Py + pz*Pz; 226 a = (a/(E + E0) - e)/E0; 226 a = (a/(E + E0) - e)/E0; 227 227 228 px = px + a*Px; 228 px = px + a*Px; 229 py = py + a*Py; 229 py = py + a*Py; 230 pz = pz + a*Pz; 230 pz = pz + a*Pz; 231 G4double p2 = px*px + py*py + pz*pz; 231 G4double p2 = px*px + py*py + pz*pz; 232 pp = std::sqrt(p2); 232 pp = std::sqrt(p2); 233 e = std::sqrt(e0*e0 + p2); 233 e = std::sqrt(e0*e0 + p2); 234 G4double ekin = e - theSecondary->GetPar 234 G4double ekin = e - theSecondary->GetParticle()->GetDefinition()->GetPDGMass()/MeV; 235 theSecondary->GetParticle()->SetMomentum 235 theSecondary->GetParticle()->SetMomentumDirection(G4ParticleMomentum(px/pp, 236 236 py/pp, 237 237 pz/pp)); 238 theSecondary->GetParticle()->SetKineticE 238 theSecondary->GetParticle()->SetKineticEnergy(ekin*MeV); 239 } 239 } 240 240 241 return &theParticleChange; 241 return &theParticleChange; 242 } 242 } 243 243 244 // Computes atomic mass in MeV (translation of 244 // Computes atomic mass in MeV (translation of GHEISHA routine ATOMAS) 245 // Not optimized: conforms closely to original 245 // Not optimized: conforms closely to original Fortran. 246 246 247 G4double G4LFission::Atomas(const G4double A, 247 G4double G4LFission::Atomas(const G4double A, const G4double Z) 248 { 248 { 249 G4double rmel = G4Electron::ElectronDefiniti 249 G4double rmel = G4Electron::ElectronDefinition()->GetPDGMass()/MeV; 250 G4double rmp = G4Proton::ProtonDefinition() 250 G4double rmp = G4Proton::ProtonDefinition()->GetPDGMass()/MeV; 251 G4double rmn = G4Neutron::NeutronDefinition 251 G4double rmn = G4Neutron::NeutronDefinition()->GetPDGMass()/MeV; 252 G4double rmd = G4Deuteron::DeuteronDefiniti 252 G4double rmd = G4Deuteron::DeuteronDefinition()->GetPDGMass()/MeV; 253 G4double rma = G4Alpha::AlphaDefinition()-> 253 G4double rma = G4Alpha::AlphaDefinition()->GetPDGMass()/MeV; 254 254 255 G4int ia = static_cast<G4int>(A + 0.5); 255 G4int ia = static_cast<G4int>(A + 0.5); 256 if (ia < 1) return 0; 256 if (ia < 1) return 0; 257 G4int iz = static_cast<G4int>(Z + 0.5); 257 G4int iz = static_cast<G4int>(Z + 0.5); 258 if (iz < 0) return 0; 258 if (iz < 0) return 0; 259 if (iz > ia) return 0; 259 if (iz > ia) return 0; 260 260 261 if (ia == 1) { 261 if (ia == 1) { 262 if (iz == 0) return rmn; //neut 262 if (iz == 0) return rmn; //neutron 263 if (iz == 1) return rmp + rmel; //Hydr 263 if (iz == 1) return rmp + rmel; //Hydrogen 264 } 264 } 265 else if (ia == 2 && iz == 1) { 265 else if (ia == 2 && iz == 1) { 266 return rmd; //Deut 266 return rmd; //Deuteron 267 } 267 } 268 else if (ia == 4 && iz == 2) { 268 else if (ia == 4 && iz == 2) { 269 return rma; //Alph 269 return rma; //Alpha 270 } 270 } 271 271 272 G4Pow* Pow=G4Pow::GetInstance(); 272 G4Pow* Pow=G4Pow::GetInstance(); 273 G4double mass = (A - Z)*rmn + Z*rmp + Z*rmel 273 G4double mass = (A - Z)*rmn + Z*rmp + Z*rmel - 15.67*A 274 + 17.23*Pow->A23(A) 274 + 17.23*Pow->A23(A) 275 + 93.15*(A/2. - Z)*(A/2. - Z 275 + 93.15*(A/2. - Z)*(A/2. - Z)/A 276 + 0.6984523*Z*Z/Pow->A13(A); 276 + 0.6984523*Z*Z/Pow->A13(A); 277 G4int ipp = (ia - iz)%2; 277 G4int ipp = (ia - iz)%2; 278 G4int izz = iz%2; 278 G4int izz = iz%2; 279 if (ipp == izz) mass = mass + (ipp + izz -1) 279 if (ipp == izz) mass = mass + (ipp + izz -1)*12.*Pow->powA(A, -0.5); 280 280 281 return mass; 281 return mass; 282 } 282 } 283 283 284 const std::pair<G4double, G4double> G4LFission 284 const std::pair<G4double, G4double> G4LFission::GetFatalEnergyCheckLevels() const 285 { 285 { 286 // max energy non-conservation is mass of he 286 // max energy non-conservation is mass of heavy nucleus 287 return std::pair<G4double, G4double>(10.0*pe 287 return std::pair<G4double, G4double>(10.0*perCent, 350.0*CLHEP::GeV); 288 } 288 } 289 289