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