- Timestamp:
- May 28, 2009, 4:26:57 PM (15 years ago)
- Location:
- trunk/source/processes/hadronic/models/de_excitation/fission/src
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/source/processes/hadronic/models/de_excitation/fission/src/G4CompetitiveFission.cc
r1007 r1055 25 25 // 26 26 // 27 // $Id: G4CompetitiveFission.cc,v 1. 9 2008/11/20 13:46:27 dennisExp $28 // GEANT4 tag $Name: geant4-09-0 2$27 // $Id: G4CompetitiveFission.cc,v 1.11 2009/03/13 18:57:17 vnivanch Exp $ 28 // GEANT4 tag $Name: geant4-09-03-beta-cand-01 $ 29 29 // 30 30 // Hadronic Process: Nuclear De-excitations 31 31 // by V. Lara (Oct 1998) 32 // 33 // J. M. Quesada (March 2009). Bugs fixed: 34 // - Full relativistic calculation (Lorentz boosts) 35 // - Fission pairing energy is included in fragment excitation energies 36 // Now Energy and momentum are conserved in fission 32 37 33 38 #include "G4CompetitiveFission.hh" … … 127 132 } 128 133 134 129 135 // Atomic Mass of Nucleus (in MeV) 130 G4double M = G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(Z,A)/MeV; 136 G4double M = G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(Z,A); 137 131 138 // Nucleus Momentum 132 139 G4LorentzVector theNucleusMomentum = theNucleus.GetMomentum(); … … 147 154 G4double FragmentsExcitationEnergy = 0.0; 148 155 G4double FragmentsKineticEnergy = 0.0; 156 157 //JMQ 04/03/09 It will be used latter to fix the bug in energy conservation 158 G4double FissionPairingEnergy=G4PairingCorrection::GetInstance()->GetFissionPairingCorrection(A,Z); 149 159 150 160 G4int Trials = 0; … … 162 172 if (A2 < 1 || Z2 < 0) 163 173 throw G4HadronicException(__FILE__, __LINE__, "G4CompetitiveFission::BreakUp: Can't define second fragment! "); 164 M2 = G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(Z2,A2)/MeV; 174 175 M2 = G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(Z2,A2); 165 176 166 177 // Check that fragment masses are less or equal than total energy 167 // if (M1 + M2 > theNucleusMomentum.mag()/MeV) 168 if (M1 + M2 > theNucleusMomentum.e()/MeV) 178 if (M1 + M2 > theNucleusMomentum.e()) 169 179 throw G4HadronicException(__FILE__, __LINE__, "G4CompetitiveFission::BreakUp: Fragments Mass > Total Energy"); 170 180 171 181 // Maximal Kinetic Energy (available energy for fragments) 172 // G4double Tmax = theNucleusMomentum.mag()/MeV - M1 - M2;173 182 G4double Tmax = M + U - M1 - M2; 174 183 … … 180 189 181 190 // Excitation Energy 182 FragmentsExcitationEnergy = Tmax - FragmentsKineticEnergy; 191 // FragmentsExcitationEnergy = Tmax - FragmentsKineticEnergy; 192 // JMQ 04/03/09 BUG FIXED: in order to fulfill energy conservation the 193 // fragments carry the fission pairing energy in form of 194 //excitation energy 195 196 FragmentsExcitationEnergy = Tmax - FragmentsKineticEnergy+FissionPairingEnergy; 197 198 } while (FragmentsExcitationEnergy < 0.0 && Trials++ < 100); 183 199 184 } while (FragmentsExcitationEnergy < 0.0 && Trials++ < 100);185 186 200 187 201 … … 190 204 191 205 192 // while (FragmentsExcitationEnergy < 0 && Trials < 100);193 194 // Fragment 1206 // while (FragmentsExcitationEnergy < 0 && Trials < 100); 207 208 // Fragment 1 195 209 G4double U1 = FragmentsExcitationEnergy * (static_cast<G4double>(A1)/static_cast<G4double>(A)); 196 210 // Fragment 2 … … 198 212 199 213 200 G4double Pmax = std::sqrt( 2 * ( ( (M1+U1)*(M2+U2) ) / 201 ( (M1+U1)+(M2+U2) ) ) * FragmentsKineticEnergy); 202 203 G4ParticleMomentum momentum1 = IsotropicVector( Pmax ); 204 G4ParticleMomentum momentum2( -momentum1 ); 214 //JMQ 04/03/09 Full relativistic calculation is performed 215 // 216 G4double Fragment1KineticEnergy=(FragmentsKineticEnergy*(FragmentsKineticEnergy+2*(M2+U2)))/(2*(M1+U1+M2+U2+FragmentsKineticEnergy)); 217 G4ParticleMomentum Momentum1(IsotropicVector(std::sqrt(Fragment1KineticEnergy*(Fragment1KineticEnergy+2*(M1+U1))))); 218 G4ParticleMomentum Momentum2(-Momentum1); 219 G4LorentzVector FourMomentum1(Momentum1,std::sqrt(Momentum1.mag2()+(M1+U1)*(M1+U1))); 220 G4LorentzVector FourMomentum2(Momentum2,std::sqrt(Momentum2.mag2()+(M2+U2)*(M2+U2))); 221 222 //JMQ 04/03/09 now we do Lorentz boosts (instead of Galileo boosts) 223 FourMomentum1.boost(theNucleusMomentum.boostVector()); 224 FourMomentum2.boost(theNucleusMomentum.boostVector()); 225 226 //////////JMQ 04/03: Old version calculation is commented 227 // There was vioation of energy momentum conservation 228 229 // G4double Pmax = std::sqrt( 2 * ( ( (M1+U1)*(M2+U2) ) / 230 // ( (M1+U1)+(M2+U2) ) ) * FragmentsKineticEnergy); 231 232 //G4ParticleMomentum momentum1 = IsotropicVector( Pmax ); 233 // G4ParticleMomentum momentum2( -momentum1 ); 205 234 206 235 // Perform a Galileo boost for fragments 207 momentum1 += (theNucleusMomentum.boostVector() * (M1+U1));208 momentum2 += (theNucleusMomentum.boostVector() * (M2+U2));236 // momentum1 += (theNucleusMomentum.boostVector() * (M1+U1)); 237 // momentum2 += (theNucleusMomentum.boostVector() * (M2+U2)); 209 238 210 239 211 240 // Create 4-momentum for first fragment 212 241 // Warning!! Energy conservation is broken 213 G4LorentzVector FourMomentum1( momentum1 , std::sqrt(momentum1.mag2() + (M1+U1)*(M1+U1))); 242 //JMQ 04/03/09 ...NOT ANY MORE!! BUGS FIXED: Energy and momentum are NOW conserved 243 // G4LorentzVector FourMomentum1( momentum1 , std::sqrt(momentum1.mag2() + (M1+U1)*(M1+U1))); 214 244 215 245 // Create 4-momentum for second fragment 216 246 // Warning!! Energy conservation is broken 217 G4LorentzVector FourMomentum2( momentum2 , std::sqrt(momentum2.mag2() + (M2+U2)*(M2+U2))); 247 //JMQ 04/03/09 ...NOT ANY MORE!! BUGS FIXED: Energy and momentum are NOW conserved 248 // G4LorentzVector FourMomentum2( momentum2 , std::sqrt(momentum2.mag2() + (M2+U2)*(M2+U2))); 249 250 ////////// 218 251 219 252 // Create Fragments … … 227 260 Fragment2->SetCreatorModel(G4String("G4CompetitiveFission")); 228 261 #endif 229 // Create Fragment Vector262 // Create Fragment Vector 230 263 G4FragmentVector * theResult = new G4FragmentVector; 231 264 -
trunk/source/processes/hadronic/models/de_excitation/fission/src/G4FissionProbability.cc
r1007 r1055 25 25 // 26 26 // 27 // $Id: G4FissionProbability.cc,v 1. 8 2007/02/12 09:39:58 ahowardExp $28 // GEANT4 tag $Name: geant4-09-0 2$27 // $Id: G4FissionProbability.cc,v 1.9 2009/02/15 17:03:25 vnivanch Exp $ 28 // GEANT4 tag $Name: geant4-09-03-beta-cand-01 $ 29 29 // 30 30 // Hadronic Process: Nuclear De-excitations 31 31 // by V. Lara (Oct 1998) 32 32 // 33 // 34 // J.M.Quesada (14 february 2009) bug fixed in fission width: missing parenthesis in the denominator 33 35 34 36 … … 95 97 G4double Exp1 = 0.0; 96 98 if (SystemEntropy <= 160.0) Exp1 = std::exp(-SystemEntropy); 97 // @@@@@@@@@@@@@@@@@ hpw changed max to min - cannot notify vicente now since cern mail gave up on me...99 // @@@@@@@@@@@@@@@@@ hpw changed max to min - cannot notify vicente now 98 100 G4double Exp2 = std::exp( std::min(700.0,Cf-SystemEntropy) ); 99 101 102 // JMQ 14/02/09 BUG fixed in fission probability (missing parenthesis at denominator) 100 103 //AH fix from Vincente: G4double probability = (Exp1 + (1.0-Cf)*Exp2) / 4.0*pi*afission; 101 G4double probability = (Exp1 + (Cf-1.0)*Exp2) / 4.0*pi*afission; 104 // G4double probability = (Exp1 + (Cf-1.0)*Exp2) / 4.0*pi*afission; 105 G4double probability = (Exp1 + (Cf-1.0)*Exp2) / (4.0*pi*afission); 106 107 102 108 103 109 return probability;
Note: See TracChangeset
for help on using the changeset viewer.