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