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 // $Id: G4CompetitiveFission.cc,v 1.14 2010-11-17 20:22:46 vnivanch Exp $ >> 28 // GEANT4 tag $Name: geant4-09-04-patch-01 $ >> 29 // 27 // Hadronic Process: Nuclear De-excitations 30 // Hadronic Process: Nuclear De-excitations 28 // by V. Lara (Oct 1998) 31 // by V. Lara (Oct 1998) 29 // 32 // 30 // J. M. Quesada (March 2009). Bugs fixed: 33 // J. M. Quesada (March 2009). Bugs fixed: 31 // - Full relativistic calculation (L 34 // - Full relativistic calculation (Lorentz boosts) 32 // - Fission pairing energy is includ 35 // - Fission pairing energy is included in fragment excitation energies 33 // Now Energy and momentum are conserved in fi 36 // Now Energy and momentum are conserved in fission 34 37 35 #include "G4CompetitiveFission.hh" 38 #include "G4CompetitiveFission.hh" 36 #include "G4PairingCorrection.hh" 39 #include "G4PairingCorrection.hh" 37 #include "G4ParticleMomentum.hh" 40 #include "G4ParticleMomentum.hh" 38 #include "G4NuclearLevelData.hh" << 39 #include "G4VFissionBarrier.hh" << 40 #include "G4FissionBarrier.hh" << 41 #include "G4FissionProbability.hh" << 42 #include "G4VLevelDensityParameter.hh" << 43 #include "G4FissionLevelDensityParameter.hh" << 44 #include "G4Pow.hh" 41 #include "G4Pow.hh" 45 #include "Randomize.hh" << 46 #include "G4RandomDirection.hh" << 47 #include "G4PhysicalConstants.hh" << 48 #include "G4PhysicsModelCatalog.hh" << 49 << 50 G4CompetitiveFission::G4CompetitiveFission() : << 51 { << 52 theFissionBarrierPtr = new G4FissionBarrier; << 53 theFissionProbabilityPtr = new G4FissionProb << 54 theLevelDensityPtr = new G4FissionLevelDensi << 55 pairingCorrection = G4NuclearLevelData::GetI << 56 theSecID = G4PhysicsModelCatalog::GetModelID << 57 } << 58 42 59 G4CompetitiveFission::~G4CompetitiveFission() << 43 G4CompetitiveFission::G4CompetitiveFission() : G4VEvaporationChannel("fission") 60 { 44 { 61 if (myOwnFissionBarrier) delete theFissionBa << 45 theFissionBarrierPtr = new G4FissionBarrier; 62 if (myOwnFissionProbability) delete theFissi << 46 MyOwnFissionBarrier = true; 63 if (myOwnLevelDensity) delete theLevelDensit << 47 >> 48 theFissionProbabilityPtr = new G4FissionProbability; >> 49 MyOwnFissionProbability = true; >> 50 >> 51 theLevelDensityPtr = new G4FissionLevelDensityParameter; >> 52 MyOwnLevelDensity = true; >> 53 >> 54 MaximalKineticEnergy = -1000.0*MeV; >> 55 FissionBarrier = 0.0; >> 56 FissionProbability = 0.0; >> 57 LevelDensityParameter = 0.0; 64 } 58 } 65 59 66 void G4CompetitiveFission::Initialise() << 60 G4CompetitiveFission::~G4CompetitiveFission() 67 { 61 { 68 if (!isInitialised) { << 62 if (MyOwnFissionBarrier) delete theFissionBarrierPtr; 69 isInitialised = true; << 63 70 G4VEvaporationChannel::Initialise(); << 64 if (MyOwnFissionProbability) delete theFissionProbabilityPtr; 71 if (OPTxs == 1) { fFactor = 0.5; } << 65 72 } << 66 if (MyOwnLevelDensity) delete theLevelDensityPtr; 73 } 67 } 74 68 75 G4double G4CompetitiveFission::GetEmissionProb << 69 void G4CompetitiveFission::Initialize(const G4Fragment & fragment) 76 { 70 { 77 if (!isInitialised) { Initialise(); } << 71 G4int anA = fragment.GetA_asInt(); 78 G4int Z = fragment->GetZ_asInt(); << 72 G4int aZ = fragment.GetZ_asInt(); 79 G4int A = fragment->GetA_asInt(); << 73 G4double ExEnergy = fragment.GetExcitationEnergy() - 80 fissionProbability = 0.0; << 74 G4PairingCorrection::GetInstance()->GetFissionPairingCorrection(anA,aZ); 81 // Saddle point excitation energy ---> A = 6 << 82 if (A >= 65 && Z > 16) { << 83 G4double exEnergy = fragment->GetExcitatio << 84 pairingCorrection->GetFissionPairingCorr << 85 75 86 if (exEnergy > 0.0) { << 76 87 fissionBarrier = theFissionBarrierPtr->F << 77 // Saddle point excitation energy ---> A = 65 88 maxKineticEnergy = exEnergy - fissionBar << 78 // Fission is excluded for A < 65 89 fissionProbability = << 79 if (anA >= 65 && ExEnergy > 0.0) { 90 theFissionProbabilityPtr->EmissionProbabilit << 80 FissionBarrier = theFissionBarrierPtr->FissionBarrier(anA,aZ,ExEnergy); 91 maxKineticEnergy); << 81 MaximalKineticEnergy = ExEnergy - FissionBarrier; >> 82 LevelDensityParameter = >> 83 theLevelDensityPtr->LevelDensityParameter(anA,aZ,ExEnergy); >> 84 FissionProbability = >> 85 theFissionProbabilityPtr->EmissionProbability(fragment,MaximalKineticEnergy); 92 } 86 } >> 87 else { >> 88 MaximalKineticEnergy = -1000.0*MeV; >> 89 LevelDensityParameter = 0.0; >> 90 FissionProbability = 0.0; 93 } 91 } 94 return fissionProbability*fFactor; << 95 } 92 } 96 93 97 G4Fragment* G4CompetitiveFission::EmittedFragm << 94 G4FragmentVector * G4CompetitiveFission::BreakUp(const G4Fragment & theNucleus) 98 { 95 { 99 G4Fragment * Fragment1 = nullptr; << 100 // Nucleus data 96 // Nucleus data 101 // Atomic number of nucleus 97 // Atomic number of nucleus 102 G4int A = theNucleus->GetA_asInt(); << 98 G4int A = theNucleus.GetA_asInt(); 103 // Charge of nucleus 99 // Charge of nucleus 104 G4int Z = theNucleus->GetZ_asInt(); << 100 G4int Z = theNucleus.GetZ_asInt(); 105 // Excitation energy (in MeV) 101 // Excitation energy (in MeV) 106 G4double U = theNucleus->GetExcitationEnergy << 102 G4double U = theNucleus.GetExcitationEnergy() - 107 G4double pcorr = pairingCorrection->GetFissi << 103 G4PairingCorrection::GetInstance()->GetFissionPairingCorrection(A,Z); 108 if (U <= pcorr) { return Fragment1; } << 104 // Check that U > 0 >> 105 if (U <= 0.0) { >> 106 G4FragmentVector * theResult = new G4FragmentVector; >> 107 theResult->push_back(new G4Fragment(theNucleus)); >> 108 return theResult; >> 109 } 109 110 110 // Atomic Mass of Nucleus (in MeV) 111 // Atomic Mass of Nucleus (in MeV) 111 G4double M = theNucleus->GetGroundStateMass( << 112 G4double M = theNucleus.GetGroundStateMass(); 112 113 113 // Nucleus Momentum 114 // Nucleus Momentum 114 G4LorentzVector theNucleusMomentum = theNucl << 115 G4LorentzVector theNucleusMomentum = theNucleus.GetMomentum(); 115 116 116 // Calculate fission parameters 117 // Calculate fission parameters 117 theParam.DefineParameters(A, Z, U-pcorr, fis << 118 G4FissionParameters theParameters(A,Z,U,FissionBarrier); 118 119 119 // First fragment 120 // First fragment 120 G4int A1 = 0; 121 G4int A1 = 0; 121 G4int Z1 = 0; 122 G4int Z1 = 0; 122 G4double M1 = 0.0; 123 G4double M1 = 0.0; 123 124 124 // Second fragment 125 // Second fragment 125 G4int A2 = 0; 126 G4int A2 = 0; 126 G4int Z2 = 0; 127 G4int Z2 = 0; 127 G4double M2 = 0.0; 128 G4double M2 = 0.0; 128 129 129 G4double FragmentsExcitationEnergy = 0.0; 130 G4double FragmentsExcitationEnergy = 0.0; 130 G4double FragmentsKineticEnergy = 0.0; 131 G4double FragmentsKineticEnergy = 0.0; 131 132 >> 133 //JMQ 04/03/09 It will be used latter to fix the bug in energy conservation >> 134 G4double FissionPairingEnergy= >> 135 G4PairingCorrection::GetInstance()->GetFissionPairingCorrection(A,Z); >> 136 132 G4int Trials = 0; 137 G4int Trials = 0; 133 do { 138 do { 134 139 135 // First fragment 140 // First fragment 136 A1 = FissionAtomicNumber(A); << 141 A1 = FissionAtomicNumber(A,theParameters); 137 Z1 = FissionCharge(A, Z, A1); << 142 Z1 = FissionCharge(A,Z,A1); 138 M1 = G4NucleiProperties::GetNuclearMass(A1 << 143 M1 = G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(Z1,A1); 139 144 140 // Second Fragment 145 // Second Fragment 141 A2 = A - A1; 146 A2 = A - A1; 142 Z2 = Z - Z1; 147 Z2 = Z - Z1; 143 if (A2 < 1 || Z2 < 0 || Z2 > A2) { << 148 if (A2 < 1 || Z2 < 0) { 144 FragmentsExcitationEnergy = -1.0; << 149 throw G4HadronicException(__FILE__, __LINE__, 145 continue; << 150 "G4CompetitiveFission::BreakUp: Can't define second fragment! "); 146 } 151 } 147 M2 = G4NucleiProperties::GetNuclearMass(A2 << 152 M2 = G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(Z2,A2); 148 // Maximal Kinetic Energy (available energ << 149 G4double Tmax = M + U - M1 - M2 - pcorr; << 150 153 151 // Check that fragment masses are less or 154 // Check that fragment masses are less or equal than total energy 152 if (Tmax < 0.0) { << 155 if (M1 + M2 > theNucleusMomentum.e()) { 153 FragmentsExcitationEnergy = -1.0; << 156 throw G4HadronicException(__FILE__, __LINE__, 154 continue; << 157 "G4CompetitiveFission::BreakUp: Fragments Mass > Total Energy"); 155 } 158 } >> 159 // Maximal Kinetic Energy (available energy for fragments) >> 160 G4double Tmax = M + U - M1 - M2; 156 161 157 FragmentsKineticEnergy = FissionKineticEne 162 FragmentsKineticEnergy = FissionKineticEnergy( A , Z, 158 A1, Z1, 163 A1, Z1, 159 A2, Z2, 164 A2, Z2, 160 U , Tmax); << 165 U , Tmax, >> 166 theParameters); 161 167 162 // Excitation Energy 168 // Excitation Energy 163 // FragmentsExcitationEnergy = Tmax - Frag << 169 // FragmentsExcitationEnergy = Tmax - FragmentsKineticEnergy; 164 // JMQ 04/03/09 BUG FIXED: in order to ful 170 // JMQ 04/03/09 BUG FIXED: in order to fulfill energy conservation the 165 // fragments carry the fission pairing ene 171 // fragments carry the fission pairing energy in form of 166 // excitation energy << 172 //excitation energy 167 173 168 FragmentsExcitationEnergy = 174 FragmentsExcitationEnergy = 169 Tmax - FragmentsKineticEnergy + pcorr; << 175 Tmax - FragmentsKineticEnergy+FissionPairingEnergy; 170 176 171 // Loop checking, 05-Aug-2015, Vladimir Iv << 177 } while (FragmentsExcitationEnergy < 0.0 && Trials++ < 100); 172 } while (FragmentsExcitationEnergy < 0.0 && << 173 178 174 if (FragmentsExcitationEnergy <= 0.0) { 179 if (FragmentsExcitationEnergy <= 0.0) { 175 throw G4HadronicException(__FILE__, __LINE 180 throw G4HadronicException(__FILE__, __LINE__, 176 "G4CompetitiveFission::BreakItUp: Excita 181 "G4CompetitiveFission::BreakItUp: Excitation energy for fragments < 0.0!"); 177 } 182 } 178 183 >> 184 // while (FragmentsExcitationEnergy < 0 && Trials < 100); >> 185 179 // Fragment 1 186 // Fragment 1 180 M1 += FragmentsExcitationEnergy * A1/static_ << 187 G4double U1 = FragmentsExcitationEnergy * A1/static_cast<G4double>(A); 181 // Fragment 2 << 188 // Fragment 2 182 M2 += FragmentsExcitationEnergy * A2/static_ << 189 G4double U2 = FragmentsExcitationEnergy * A2/static_cast<G4double>(A); 183 // primary << 190 184 M += U; << 191 //JMQ 04/03/09 Full relativistic calculation is performed 185 << 192 // 186 G4double etot1 = ((M - M2)*(M + M2) + M1*M1) << 193 G4double Fragment1KineticEnergy= 187 G4ParticleMomentum Momentum1 = << 194 (FragmentsKineticEnergy*(FragmentsKineticEnergy+2*(M2+U2))) 188 std::sqrt((etot1 - M1)*(etot1+M1))*G4Rando << 195 /(2*(M1+U1+M2+U2+FragmentsKineticEnergy)); 189 G4LorentzVector FourMomentum1(Momentum1, eto << 196 G4ParticleMomentum Momentum1(IsotropicVector(std::sqrt(Fragment1KineticEnergy*(Fragment1KineticEnergy+2*(M1+U1))))); >> 197 G4ParticleMomentum Momentum2(-Momentum1); >> 198 G4LorentzVector FourMomentum1(Momentum1,std::sqrt(Momentum1.mag2()+(M1+U1)*(M1+U1))); >> 199 G4LorentzVector FourMomentum2(Momentum2,std::sqrt(Momentum2.mag2()+(M2+U2)*(M2+U2))); >> 200 >> 201 //JMQ 04/03/09 now we do Lorentz boosts (instead of Galileo boosts) 190 FourMomentum1.boost(theNucleusMomentum.boost 202 FourMomentum1.boost(theNucleusMomentum.boostVector()); >> 203 FourMomentum2.boost(theNucleusMomentum.boostVector()); 191 204 >> 205 //////////JMQ 04/03: Old version calculation is commented >> 206 // There was vioation of energy momentum conservation >> 207 >> 208 // G4double Pmax = std::sqrt( 2 * ( ( (M1+U1)*(M2+U2) ) / >> 209 // ( (M1+U1)+(M2+U2) ) ) * FragmentsKineticEnergy); >> 210 >> 211 //G4ParticleMomentum momentum1 = IsotropicVector( Pmax ); >> 212 // G4ParticleMomentum momentum2( -momentum1 ); >> 213 >> 214 // Perform a Galileo boost for fragments >> 215 // momentum1 += (theNucleusMomentum.boostVector() * (M1+U1)); >> 216 // momentum2 += (theNucleusMomentum.boostVector() * (M2+U2)); >> 217 >> 218 >> 219 // Create 4-momentum for first fragment >> 220 // Warning!! Energy conservation is broken >> 221 //JMQ 04/03/09 ...NOT ANY MORE!! BUGS FIXED: Energy and momentum are NOW conserved >> 222 // G4LorentzVector FourMomentum1( momentum1 , std::sqrt(momentum1.mag2() + (M1+U1)*(M1+U1))); >> 223 >> 224 // Create 4-momentum for second fragment >> 225 // Warning!! Energy conservation is broken >> 226 //JMQ 04/03/09 ...NOT ANY MORE!! BUGS FIXED: Energy and momentum are NOW conserved >> 227 // G4LorentzVector FourMomentum2( momentum2 , std::sqrt(momentum2.mag2() + (M2+U2)*(M2+U2))); >> 228 >> 229 ////////// >> 230 192 // Create Fragments 231 // Create Fragments 193 Fragment1 = new G4Fragment( A1, Z1, FourMome << 232 G4Fragment * Fragment1 = new G4Fragment( A1, Z1, FourMomentum1); 194 if (Fragment1 != nullptr) { Fragment1->SetCr << 233 G4Fragment * Fragment2 = new G4Fragment( A2, Z2, FourMomentum2); 195 theNucleusMomentum -= FourMomentum1; << 234 196 theNucleus->SetZandA_asInt(Z2, A2); << 235 // Create Fragment Vector 197 theNucleus->SetMomentum(theNucleusMomentum); << 236 G4FragmentVector * theResult = new G4FragmentVector; 198 theNucleus->SetCreatorModelID(theSecID); << 237 199 return Fragment1; << 238 theResult->push_back(Fragment1); >> 239 theResult->push_back(Fragment2); >> 240 >> 241 #ifdef debug >> 242 CheckConservation(theNucleus,theResult); >> 243 #endif >> 244 >> 245 return theResult; 200 } 246 } 201 247 202 G4int 248 G4int 203 G4CompetitiveFission::FissionAtomicNumber(G4in << 249 G4CompetitiveFission::FissionAtomicNumber(G4int A, >> 250 const G4FissionParameters & theParam) 204 // Calculates the atomic number of a fission 251 // Calculates the atomic number of a fission product 205 { 252 { 206 253 207 // For Simplicity reading code 254 // For Simplicity reading code 208 G4int A1 = theParam.GetA1(); << 255 const G4double A1 = theParam.GetA1(); 209 G4int A2 = theParam.GetA2(); << 256 const G4double A2 = theParam.GetA2(); 210 G4double As = theParam.GetAs(); << 257 const G4double As = theParam.GetAs(); 211 G4double Sigma2 = theParam.GetSigma2(); << 258 // const G4double Sigma1 = theParam.GetSigma1(); 212 G4double SigmaS = theParam.GetSigmaS(); << 259 const G4double Sigma2 = theParam.GetSigma2(); 213 G4double w = theParam.GetW(); << 260 const G4double SigmaS = theParam.GetSigmaS(); >> 261 const G4double w = theParam.GetW(); 214 262 >> 263 // G4double FasymAsym = 2.0*std::exp(-((A2-As)*(A2-As))/(2.0*Sigma2*Sigma2)) + >> 264 // std::exp(-((A1-As)*(A1-As))/(2.0*Sigma1*Sigma1)); >> 265 >> 266 // G4double FsymA1A2 = std::exp(-((As-(A1+A2))*(As-(A1+A2)))/(2.0*SigmaS*SigmaS)); >> 267 215 G4double C2A = A2 + 3.72*Sigma2; 268 G4double C2A = A2 + 3.72*Sigma2; 216 G4double C2S = As + 3.72*SigmaS; 269 G4double C2S = As + 3.72*SigmaS; 217 270 218 G4double C2 = 0.0; 271 G4double C2 = 0.0; 219 if (w > 1000.0 ) { C2 = C2S; } << 272 if (w > 1000.0 ) C2 = C2S; 220 else if (w < 0.001) { C2 = C2A; } << 273 else if (w < 0.001) C2 = C2A; 221 else { C2 = std::max(C2A,C2S << 274 else C2 = std::max(C2A,C2S); 222 275 223 G4double C1 = A-C2; 276 G4double C1 = A-C2; 224 if (C1 < 30.0) { 277 if (C1 < 30.0) { 225 C2 = A-30.0; 278 C2 = A-30.0; 226 C1 = 30.0; 279 C1 = 30.0; 227 } 280 } 228 281 229 G4double Am1 = (As + A1)*0.5; << 282 G4double Am1 = (As + A1)/2.0; 230 G4double Am2 = (A1 + A2)*0.5; << 283 G4double Am2 = (A1 + A2)/2.0; 231 284 232 // Get Mass distributions as sum of symmetri 285 // Get Mass distributions as sum of symmetric and asymmetric Gasussians 233 G4double Mass1 = MassDistribution(As,A); << 286 G4double Mass1 = MassDistribution(As,A,theParam); 234 G4double Mass2 = MassDistribution(Am1,A); << 287 G4double Mass2 = MassDistribution(Am1,A,theParam); 235 G4double Mass3 = MassDistribution(G4double(A << 288 G4double Mass3 = MassDistribution(A1,A,theParam); 236 G4double Mass4 = MassDistribution(Am2,A); << 289 G4double Mass4 = MassDistribution(Am2,A,theParam); 237 G4double Mass5 = MassDistribution(G4double(A << 290 G4double Mass5 = MassDistribution(A2,A,theParam); 238 // get maximal value among Mass1,...,Mass5 291 // get maximal value among Mass1,...,Mass5 239 G4double MassMax = Mass1; 292 G4double MassMax = Mass1; 240 if (Mass2 > MassMax) { MassMax = Mass2; } << 293 if (Mass2 > MassMax) MassMax = Mass2; 241 if (Mass3 > MassMax) { MassMax = Mass3; } << 294 if (Mass3 > MassMax) MassMax = Mass3; 242 if (Mass4 > MassMax) { MassMax = Mass4; } << 295 if (Mass4 > MassMax) MassMax = Mass4; 243 if (Mass5 > MassMax) { MassMax = Mass5; } << 296 if (Mass5 > MassMax) MassMax = Mass5; 244 297 245 // Sample a fragment mass number, which lies 298 // Sample a fragment mass number, which lies between C1 and C2 246 G4double xm; << 299 G4double m; 247 G4double Pm; 300 G4double Pm; 248 do { 301 do { 249 xm = C1+G4UniformRand()*(C2-C1); << 302 m = C1+G4UniformRand()*(C2-C1); 250 Pm = MassDistribution(xm,A); << 303 Pm = MassDistribution(m,A,theParam); 251 // Loop checking, 05-Aug-2015, Vladimir Iv << 252 } while (MassMax*G4UniformRand() > Pm); 304 } while (MassMax*G4UniformRand() > Pm); 253 G4int ires = G4lrint(xm); << 254 305 255 return ires; << 306 return static_cast<G4int>(m+0.5); 256 } 307 } 257 308 258 G4double 309 G4double 259 G4CompetitiveFission::MassDistribution(G4doubl << 310 G4CompetitiveFission::MassDistribution(G4double x, G4double A, >> 311 const G4FissionParameters & theParam) 260 // This method gives mass distribution F(x) 312 // This method gives mass distribution F(x) = F_{asym}(x)+w*F_{sym}(x) 261 // which consist of symmetric and asymmetric 313 // which consist of symmetric and asymmetric sum of gaussians components. 262 { 314 { 263 G4double y0 = (x-theParam.GetAs())/theParam. << 315 G4double Xsym = std::exp(-0.5*(x-theParam.GetAs())*(x-theParam.GetAs())/ 264 G4double Xsym = LocalExp(y0); << 316 (theParam.GetSigmaS()*theParam.GetSigmaS())); 265 317 266 G4double y1 = (x - theParam.GetA1())/thePara << 318 G4double Xasym = std::exp(-0.5*(x-theParam.GetA2())*(x-theParam.GetA2())/ 267 G4double y2 = (x - theParam.GetA2())/thePara << 319 (theParam.GetSigma2()*theParam.GetSigma2())) + 268 G4double z1 = (x - A + theParam.GetA1())/the << 320 std::exp(-0.5*(x-(A-theParam.GetA2()))*(x-(A-theParam.GetA2()))/ 269 G4double z2 = (x - A + theParam.GetA2())/the << 321 (theParam.GetSigma2()*theParam.GetSigma2())) + 270 G4double Xasym = LocalExp(y1) + LocalExp(y2) << 322 0.5*std::exp(-0.5*(x-theParam.GetA1())*(x-theParam.GetA1())/ 271 + 0.5*(LocalExp(z1) + LocalExp(z2)); << 323 (theParam.GetSigma1()*theParam.GetSigma1())) + 272 << 324 0.5*std::exp(-0.5*(x-(A-theParam.GetA1()))*(x-(A-theParam.GetA1()))/ 273 G4double res; << 325 (theParam.GetSigma1()*theParam.GetSigma1())); 274 G4double w = theParam.GetW(); << 326 275 if (w > 1000) { res = Xsym; } << 327 if (theParam.GetW() > 1000) return Xsym; 276 else if (w < 0.001) { res = Xasym; } << 328 else if (theParam.GetW() < 0.001) return Xasym; 277 else { res = w*Xsym+Xasym; } << 329 else return theParam.GetW()*Xsym+Xasym; 278 return res; << 279 } 330 } 280 331 281 G4int G4CompetitiveFission::FissionCharge(G4in << 332 G4int G4CompetitiveFission::FissionCharge(G4double A, G4double Z, >> 333 G4double Af) 282 // Calculates the charge of a fission produc 334 // Calculates the charge of a fission product for a given atomic number Af 283 { 335 { 284 static const G4double sigma = 0.6; << 336 const G4double sigma = 0.6; 285 G4double DeltaZ = 0.0; 337 G4double DeltaZ = 0.0; 286 if (Af >= 134.0) { DeltaZ = -0.45; << 338 if (Af >= 134.0) DeltaZ = -0.45; // 134 <= Af 287 else if (Af <= (A-134.0)) { DeltaZ = 0.45; } << 339 else if (Af <= (A-134.0)) DeltaZ = 0.45; // Af <= (A-134) 288 else { DeltaZ = -0.45*( << 340 else DeltaZ = -0.45*(Af-(A/2.0))/(134.0-(A/2.0)); // (A-134) < Af < 134 289 341 290 G4double Zmean = (Af/A)*Z + DeltaZ; 342 G4double Zmean = (Af/A)*Z + DeltaZ; 291 343 292 G4double theZ; 344 G4double theZ; 293 do { 345 do { 294 theZ = G4RandGauss::shoot(Zmean,sigma); 346 theZ = G4RandGauss::shoot(Zmean,sigma); 295 // Loop checking, 05-Aug-2015, Vladimir Iv << 296 } while (theZ < 1.0 || theZ > (Z-1.0) || th 347 } while (theZ < 1.0 || theZ > (Z-1.0) || theZ > Af); 297 << 348 // return static_cast<G4int>(theZ+0.5); 298 return G4lrint(theZ); << 349 return static_cast<G4int>(theZ+0.5); 299 } 350 } 300 351 301 G4double 352 G4double 302 G4CompetitiveFission::FissionKineticEnergy(G4i 353 G4CompetitiveFission::FissionKineticEnergy(G4int A, G4int Z, 303 G4int Af1, G4int /*Zf1*/, << 354 G4double Af1, G4double /*Zf1*/, 304 G4int Af2, G4int /*Zf2*/, << 355 G4double Af2, G4double /*Zf2*/, 305 G4double /*U*/, G4double Tmax) << 356 G4double /*U*/, G4double Tmax, >> 357 const G4FissionParameters & theParam) 306 // Gives the kinetic energy of fission produ 358 // Gives the kinetic energy of fission products 307 { 359 { 308 // Find maximal value of A for fragments 360 // Find maximal value of A for fragments 309 G4int AfMax = std::max(Af1,Af2); << 361 G4double AfMax = std::max(Af1,Af2); >> 362 if (AfMax < (A/2.0)) AfMax = A - AfMax; 310 363 311 // Weights for symmetric and asymmetric comp 364 // Weights for symmetric and asymmetric components 312 G4double Pas = 0.0; << 365 G4double Pas; 313 if (theParam.GetW() <= 1000) { << 366 if (theParam.GetW() > 1000) Pas = 0.0; 314 G4double x1 = (AfMax-theParam.GetA1())/the << 367 else { 315 G4double x2 = (AfMax-theParam.GetA2())/the << 368 G4double P1 = 0.5*std::exp(-0.5*(AfMax-theParam.GetA1())*(AfMax-theParam.GetA1())/ 316 Pas = 0.5*LocalExp(x1) + LocalExp(x2); << 369 (theParam.GetSigma1()*theParam.GetSigma1())); >> 370 >> 371 G4double P2 = std::exp(-0.5*(AfMax-theParam.GetA2())*(AfMax-theParam.GetA2())/ >> 372 (theParam.GetSigma2()*theParam.GetSigma2())); >> 373 >> 374 Pas = P1+P2; 317 } 375 } 318 376 319 G4double Ps = 0.0; << 377 G4double Ps; 320 if (theParam.GetW() >= 0.001) { << 378 if (theParam.GetW() < 0.001) Ps = 0.0; 321 G4double xs = (AfMax-theParam.GetAs())/the << 379 else { 322 Ps = theParam.GetW()*LocalExp(xs); << 380 Ps = theParam.GetW()*std::exp(-0.5*(AfMax-theParam.GetAs())*(AfMax-theParam.GetAs())/ >> 381 (theParam.GetSigmaS()*theParam.GetSigmaS())); 323 } 382 } 324 G4double Psy = (Pas + Ps > 0.0) ? Ps/(Pas+Ps << 383 G4double Psy = Ps/(Pas+Ps); 325 384 326 // Fission fractions Xsy and Xas formed in s 385 // Fission fractions Xsy and Xas formed in symmetric and asymmetric modes 327 G4double PPas = theParam.GetSigma1() + 2.0 * 386 G4double PPas = theParam.GetSigma1() + 2.0 * theParam.GetSigma2(); 328 G4double PPsy = theParam.GetW() * theParam.G 387 G4double PPsy = theParam.GetW() * theParam.GetSigmaS(); 329 G4double Xas = (PPas + PPsy > 0.0) ? PPas/(P << 388 G4double Xas = PPas / (PPas+PPsy); 330 G4double Xsy = 1.0 - Xas; << 389 G4double Xsy = PPsy / (PPas+PPsy); 331 390 332 // Average kinetic energy for symmetric and 391 // Average kinetic energy for symmetric and asymmetric components 333 G4double Eaverage = (0.1071*(Z*Z)/G4Pow::Get << 392 G4double Eaverage = 0.1071*MeV*(Z*Z)/G4Pow::GetInstance()->Z13(A) + 22.2*MeV; >> 393 334 394 335 // Compute maximal average kinetic energy of << 395 // Compute maximal average kinetic energy of fragments and Energy Dispersion (sqrt) 336 G4double TaverageAfMax; 396 G4double TaverageAfMax; 337 G4double ESigma = 10*CLHEP::MeV; << 397 G4double ESigma; 338 // Select randomly fission mode (symmetric o 398 // Select randomly fission mode (symmetric or asymmetric) 339 if (G4UniformRand() > Psy) { // Asymmetric M 399 if (G4UniformRand() > Psy) { // Asymmetric Mode 340 G4double A11 = theParam.GetA1()-0.7979*the 400 G4double A11 = theParam.GetA1()-0.7979*theParam.GetSigma1(); 341 G4double A12 = theParam.GetA1()+0.7979*the 401 G4double A12 = theParam.GetA1()+0.7979*theParam.GetSigma1(); 342 G4double A21 = theParam.GetA2()-0.7979*the 402 G4double A21 = theParam.GetA2()-0.7979*theParam.GetSigma2(); 343 G4double A22 = theParam.GetA2()+0.7979*the 403 G4double A22 = theParam.GetA2()+0.7979*theParam.GetSigma2(); 344 // scale factor 404 // scale factor 345 G4double ScaleFactor = 0.5*theParam.GetSig << 405 G4double ScaleFactor = 0.5*theParam.GetSigma1()*(AsymmetricRatio(A,A11)+AsymmetricRatio(A,A12))+ 346 (AsymmetricRatio(A,A11)+AsymmetricRatio( << 347 theParam.GetSigma2()*(AsymmetricRatio(A, 406 theParam.GetSigma2()*(AsymmetricRatio(A,A21)+AsymmetricRatio(A,A22)); 348 // Compute average kinetic energy for frag 407 // Compute average kinetic energy for fragment with AfMax 349 TaverageAfMax = (Eaverage + 12.5 * Xsy) * << 408 TaverageAfMax = (Eaverage + 12.5 * Xsy) * (PPas/ScaleFactor) * AsymmetricRatio(A,AfMax); 350 AsymmetricRatio(A,G4double(AfMax)); << 409 ESigma = 10.0*MeV; // MeV 351 410 352 } else { // Symmetric Mode 411 } else { // Symmetric Mode 353 G4double As0 = theParam.GetAs() + 0.7979*t 412 G4double As0 = theParam.GetAs() + 0.7979*theParam.GetSigmaS(); >> 413 // scale factor >> 414 G4double ScaleFactor = theParam.GetW()*theParam.GetSigmaS()*SymmetricRatio(A,As0); 354 // Compute average kinetic energy for frag 415 // Compute average kinetic energy for fragment with AfMax 355 TaverageAfMax = (Eaverage - 12.5*CLHEP::Me << 416 TaverageAfMax = (Eaverage - 12.5*MeV*Xas) * (PPsy/ScaleFactor) * SymmetricRatio(A,AfMax); 356 *SymmetricRatio(A, G4double(AfMax))/Symm << 417 ESigma = 8.0*MeV; 357 ESigma = 8.0*CLHEP::MeV; << 358 } 418 } 359 419 360 // Select randomly, in accordance with Gauss << 420 361 // fragment kinetic energy << 421 // Select randomly, in accordance with Gaussian distribution, fragment kinetic energy 362 G4double KineticEnergy; 422 G4double KineticEnergy; 363 G4int i = 0; 423 G4int i = 0; 364 do { 424 do { 365 KineticEnergy = G4RandGauss::shoot(Taverag << 425 KineticEnergy = G4RandGauss::shoot(TaverageAfMax,ESigma); 366 if (++i > 100) return Eaverage; << 426 if (i++ > 100) return Eaverage; 367 // Loop checking, 05-Aug-2015, Vladimir Iv << 368 } while (KineticEnergy < Eaverage-3.72*ESigm 427 } while (KineticEnergy < Eaverage-3.72*ESigma || 369 KineticEnergy > Eaverage+3.72*ESigma || 428 KineticEnergy > Eaverage+3.72*ESigma || 370 KineticEnergy > Tmax); 429 KineticEnergy > Tmax); 371 430 372 return KineticEnergy; 431 return KineticEnergy; 373 } 432 } 374 433 375 void G4CompetitiveFission::SetFissionBarrier(G << 434 G4double G4CompetitiveFission::AsymmetricRatio(G4int A, G4double A11) >> 435 { >> 436 const G4double B1 = 23.5; >> 437 const G4double A00 = 134.0; >> 438 return Ratio(G4double(A),A11,B1,A00); >> 439 } >> 440 >> 441 G4double G4CompetitiveFission::SymmetricRatio(G4int A, G4double A11) 376 { 442 { 377 if (myOwnFissionBarrier) delete theFissionBa << 443 const G4double B1 = 5.32; 378 theFissionBarrierPtr = aBarrier; << 444 const G4double A00 = A/2.0; 379 myOwnFissionBarrier = false; << 445 return Ratio(G4double(A),A11,B1,A00); 380 } 446 } 381 447 382 void << 448 G4double G4CompetitiveFission::Ratio(G4double A, G4double A11, 383 G4CompetitiveFission::SetEmissionStrategy(G4VE << 449 G4double B1, G4double A00) 384 { 450 { 385 if (myOwnFissionProbability) delete theFissi << 451 if (A == 0.0) { 386 theFissionProbabilityPtr = aFissionProb; << 452 throw G4HadronicException(__FILE__, __LINE__, 387 myOwnFissionProbability = false; << 453 "G4CompetitiveFission::Ratio: A == 0!"); >> 454 } >> 455 if (A11 >= A/2.0 && A11 <= (A00+10.0)) { >> 456 return 1.0-B1*((A11-A00)/A)*((A11-A00)/A); >> 457 } else { >> 458 return 1.0-B1*(10.0/A)*(10.0/A)-2.0*(10.0/A)*B1*((A11-A00-10.0)/A); >> 459 } 388 } 460 } 389 461 390 void << 462 G4ThreeVector G4CompetitiveFission::IsotropicVector(const G4double Magnitude) 391 G4CompetitiveFission::SetLevelDensityParameter << 463 // Samples a isotropic random vectorwith a magnitud given by Magnitude. 392 { << 464 // By default Magnitude = 1.0 393 if (myOwnLevelDensity) delete theLevelDensit << 465 { 394 theLevelDensityPtr = aLevelDensity; << 466 G4double CosTheta = 1.0 - 2.0*G4UniformRand(); 395 myOwnLevelDensity = false; << 467 G4double SinTheta = std::sqrt(1.0 - CosTheta*CosTheta); >> 468 G4double Phi = twopi*G4UniformRand(); >> 469 G4ThreeVector Vector(Magnitude*std::cos(Phi)*SinTheta, >> 470 Magnitude*std::sin(Phi)*SinTheta, >> 471 Magnitude*CosTheta); >> 472 return Vector; >> 473 } >> 474 >> 475 #ifdef debug >> 476 void G4CompetitiveFission::CheckConservation(const G4Fragment & theInitialState, >> 477 G4FragmentVector * Result) const >> 478 { >> 479 G4double ProductsEnergy =0; >> 480 G4ThreeVector ProductsMomentum; >> 481 G4int ProductsA = 0; >> 482 G4int ProductsZ = 0; >> 483 G4FragmentVector::iterator h; >> 484 for (h = Result->begin(); h != Result->end(); h++) { >> 485 G4LorentzVector tmp = (*h)->GetMomentum(); >> 486 ProductsEnergy += tmp.e(); >> 487 ProductsMomentum += tmp.vect(); >> 488 ProductsA += static_cast<G4int>((*h)->GetA()); >> 489 ProductsZ += static_cast<G4int>((*h)->GetZ()); >> 490 } >> 491 >> 492 if (ProductsA != theInitialState.GetA()) { >> 493 G4cout << "!!!!!!!!!! Baryonic Number Conservation Violation !!!!!!!!!!" << G4endl; >> 494 G4cout << "G4CompetitiveFission.cc: Barionic Number Conservation test for fission fragments" >> 495 << G4endl; >> 496 G4cout << "Initial A = " << theInitialState.GetA() >> 497 << " Fragments A = " << ProductsA << " Diference --> " >> 498 << theInitialState.GetA() - ProductsA << G4endl; >> 499 } >> 500 if (ProductsZ != theInitialState.GetZ()) { >> 501 G4cout << "!!!!!!!!!! Charge Conservation Violation !!!!!!!!!!" << G4endl; >> 502 G4cout << "G4CompetitiveFission.cc: Charge Conservation test for fission fragments" >> 503 << G4endl; >> 504 G4cout << "Initial Z = " << theInitialState.GetZ() >> 505 << " Fragments Z = " << ProductsZ << " Diference --> " >> 506 << theInitialState.GetZ() - ProductsZ << G4endl; >> 507 } >> 508 if (std::abs(ProductsEnergy-theInitialState.GetMomentum().e()) > 1.0*keV) { >> 509 G4cout << "!!!!!!!!!! Energy Conservation Violation !!!!!!!!!!" << G4endl; >> 510 G4cout << "G4CompetitiveFission.cc: Energy Conservation test for fission fragments" >> 511 << G4endl; >> 512 G4cout << "Initial E = " << theInitialState.GetMomentum().e()/MeV << " MeV" >> 513 << " Fragments E = " << ProductsEnergy/MeV << " MeV Diference --> " >> 514 << (theInitialState.GetMomentum().e() - ProductsEnergy)/MeV << " MeV" << G4endl; >> 515 } >> 516 if (std::abs(ProductsMomentum.x()-theInitialState.GetMomentum().x()) > 1.0*keV || >> 517 std::abs(ProductsMomentum.y()-theInitialState.GetMomentum().y()) > 1.0*keV || >> 518 std::abs(ProductsMomentum.z()-theInitialState.GetMomentum().z()) > 1.0*keV) { >> 519 G4cout << "!!!!!!!!!! Momentum Conservation Violation !!!!!!!!!!" << G4endl; >> 520 G4cout << "G4CompetitiveFission.cc: Momentum Conservation test for fission fragments" >> 521 << G4endl; >> 522 G4cout << "Initial P = " << theInitialState.GetMomentum().vect() << " MeV" >> 523 << " Fragments P = " << ProductsMomentum << " MeV Diference --> " >> 524 << theInitialState.GetMomentum().vect() - ProductsMomentum << " MeV" << G4endl; >> 525 } >> 526 return; 396 } 527 } >> 528 #endif >> 529 >> 530 >> 531 397 532 398 533