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 // Hadronic Process: Nuclear De-excitations 27 // Hadronic Process: Nuclear De-excitations 28 // by V. Lara (Oct 1998) 28 // by V. Lara (Oct 1998) 29 // 29 // 30 // J. M. Quesada (March 2009). Bugs fixed: 30 // J. M. Quesada (March 2009). Bugs fixed: 31 // - Full relativistic calculation (L 31 // - Full relativistic calculation (Lorentz boosts) 32 // - Fission pairing energy is includ 32 // - Fission pairing energy is included in fragment excitation energies 33 // Now Energy and momentum are conserved in fi 33 // Now Energy and momentum are conserved in fission 34 34 35 #include "G4CompetitiveFission.hh" 35 #include "G4CompetitiveFission.hh" 36 #include "G4PairingCorrection.hh" 36 #include "G4PairingCorrection.hh" 37 #include "G4ParticleMomentum.hh" 37 #include "G4ParticleMomentum.hh" 38 #include "G4NuclearLevelData.hh" 38 #include "G4NuclearLevelData.hh" 39 #include "G4VFissionBarrier.hh" 39 #include "G4VFissionBarrier.hh" 40 #include "G4FissionBarrier.hh" 40 #include "G4FissionBarrier.hh" 41 #include "G4FissionProbability.hh" 41 #include "G4FissionProbability.hh" 42 #include "G4VLevelDensityParameter.hh" 42 #include "G4VLevelDensityParameter.hh" 43 #include "G4FissionLevelDensityParameter.hh" 43 #include "G4FissionLevelDensityParameter.hh" 44 #include "G4Pow.hh" 44 #include "G4Pow.hh" 45 #include "Randomize.hh" 45 #include "Randomize.hh" 46 #include "G4RandomDirection.hh" 46 #include "G4RandomDirection.hh" 47 #include "G4PhysicalConstants.hh" 47 #include "G4PhysicalConstants.hh" 48 #include "G4PhysicsModelCatalog.hh" 48 #include "G4PhysicsModelCatalog.hh" 49 49 50 G4CompetitiveFission::G4CompetitiveFission() : 50 G4CompetitiveFission::G4CompetitiveFission() : G4VEvaporationChannel("fission"), theSecID(-1) 51 { 51 { 52 theFissionBarrierPtr = new G4FissionBarrier; 52 theFissionBarrierPtr = new G4FissionBarrier; >> 53 myOwnFissionBarrier = true; >> 54 53 theFissionProbabilityPtr = new G4FissionProb 55 theFissionProbabilityPtr = new G4FissionProbability; >> 56 myOwnFissionProbability = true; >> 57 54 theLevelDensityPtr = new G4FissionLevelDensi 58 theLevelDensityPtr = new G4FissionLevelDensityParameter; >> 59 myOwnLevelDensity = true; >> 60 >> 61 maxKineticEnergy = fissionBarrier = fissionProbability = 0.0; 55 pairingCorrection = G4NuclearLevelData::GetI 62 pairingCorrection = G4NuclearLevelData::GetInstance()->GetPairingCorrection(); >> 63 56 theSecID = G4PhysicsModelCatalog::GetModelID 64 theSecID = G4PhysicsModelCatalog::GetModelID("model_G4CompetitiveFission"); 57 } 65 } 58 66 59 G4CompetitiveFission::~G4CompetitiveFission() 67 G4CompetitiveFission::~G4CompetitiveFission() 60 { 68 { 61 if (myOwnFissionBarrier) delete theFissionBa 69 if (myOwnFissionBarrier) delete theFissionBarrierPtr; 62 if (myOwnFissionProbability) delete theFissi 70 if (myOwnFissionProbability) delete theFissionProbabilityPtr; 63 if (myOwnLevelDensity) delete theLevelDensit 71 if (myOwnLevelDensity) delete theLevelDensityPtr; 64 } 72 } 65 73 66 void G4CompetitiveFission::Initialise() << 67 { << 68 if (!isInitialised) { << 69 isInitialised = true; << 70 G4VEvaporationChannel::Initialise(); << 71 if (OPTxs == 1) { fFactor = 0.5; } << 72 } << 73 } << 74 << 75 G4double G4CompetitiveFission::GetEmissionProb 74 G4double G4CompetitiveFission::GetEmissionProbability(G4Fragment* fragment) 76 { 75 { 77 if (!isInitialised) { Initialise(); } << 78 G4int Z = fragment->GetZ_asInt(); 76 G4int Z = fragment->GetZ_asInt(); 79 G4int A = fragment->GetA_asInt(); 77 G4int A = fragment->GetA_asInt(); 80 fissionProbability = 0.0; 78 fissionProbability = 0.0; 81 // Saddle point excitation energy ---> A = 6 79 // Saddle point excitation energy ---> A = 65 82 if (A >= 65 && Z > 16) { 80 if (A >= 65 && Z > 16) { 83 G4double exEnergy = fragment->GetExcitatio 81 G4double exEnergy = fragment->GetExcitationEnergy() - 84 pairingCorrection->GetFissionPairingCorr 82 pairingCorrection->GetFissionPairingCorrection(A, Z); 85 83 86 if (exEnergy > 0.0) { 84 if (exEnergy > 0.0) { 87 fissionBarrier = theFissionBarrierPtr->F 85 fissionBarrier = theFissionBarrierPtr->FissionBarrier(A, Z, exEnergy); 88 maxKineticEnergy = exEnergy - fissionBar 86 maxKineticEnergy = exEnergy - fissionBarrier; 89 fissionProbability = 87 fissionProbability = 90 theFissionProbabilityPtr->EmissionProbabilit 88 theFissionProbabilityPtr->EmissionProbability(*fragment, 91 maxKineticEnergy); 89 maxKineticEnergy); 92 } 90 } 93 } 91 } 94 return fissionProbability*fFactor; << 92 return fissionProbability; 95 } 93 } 96 94 97 G4Fragment* G4CompetitiveFission::EmittedFragm 95 G4Fragment* G4CompetitiveFission::EmittedFragment(G4Fragment* theNucleus) 98 { 96 { 99 G4Fragment * Fragment1 = nullptr; 97 G4Fragment * Fragment1 = nullptr; 100 // Nucleus data 98 // Nucleus data 101 // Atomic number of nucleus 99 // Atomic number of nucleus 102 G4int A = theNucleus->GetA_asInt(); 100 G4int A = theNucleus->GetA_asInt(); 103 // Charge of nucleus 101 // Charge of nucleus 104 G4int Z = theNucleus->GetZ_asInt(); 102 G4int Z = theNucleus->GetZ_asInt(); 105 // Excitation energy (in MeV) 103 // Excitation energy (in MeV) 106 G4double U = theNucleus->GetExcitationEnergy 104 G4double U = theNucleus->GetExcitationEnergy(); 107 G4double pcorr = pairingCorrection->GetFissi 105 G4double pcorr = pairingCorrection->GetFissionPairingCorrection(A,Z); 108 if (U <= pcorr) { return Fragment1; } 106 if (U <= pcorr) { return Fragment1; } 109 107 110 // Atomic Mass of Nucleus (in MeV) 108 // Atomic Mass of Nucleus (in MeV) 111 G4double M = theNucleus->GetGroundStateMass( 109 G4double M = theNucleus->GetGroundStateMass(); 112 110 113 // Nucleus Momentum 111 // Nucleus Momentum 114 G4LorentzVector theNucleusMomentum = theNucl 112 G4LorentzVector theNucleusMomentum = theNucleus->GetMomentum(); 115 113 116 // Calculate fission parameters 114 // Calculate fission parameters 117 theParam.DefineParameters(A, Z, U-pcorr, fis 115 theParam.DefineParameters(A, Z, U-pcorr, fissionBarrier); 118 116 119 // First fragment 117 // First fragment 120 G4int A1 = 0; 118 G4int A1 = 0; 121 G4int Z1 = 0; 119 G4int Z1 = 0; 122 G4double M1 = 0.0; 120 G4double M1 = 0.0; 123 121 124 // Second fragment 122 // Second fragment 125 G4int A2 = 0; 123 G4int A2 = 0; 126 G4int Z2 = 0; 124 G4int Z2 = 0; 127 G4double M2 = 0.0; 125 G4double M2 = 0.0; 128 126 129 G4double FragmentsExcitationEnergy = 0.0; 127 G4double FragmentsExcitationEnergy = 0.0; 130 G4double FragmentsKineticEnergy = 0.0; 128 G4double FragmentsKineticEnergy = 0.0; 131 129 132 G4int Trials = 0; 130 G4int Trials = 0; 133 do { 131 do { 134 132 135 // First fragment 133 // First fragment 136 A1 = FissionAtomicNumber(A); 134 A1 = FissionAtomicNumber(A); 137 Z1 = FissionCharge(A, Z, A1); 135 Z1 = FissionCharge(A, Z, A1); 138 M1 = G4NucleiProperties::GetNuclearMass(A1 136 M1 = G4NucleiProperties::GetNuclearMass(A1, Z1); 139 137 140 // Second Fragment 138 // Second Fragment 141 A2 = A - A1; 139 A2 = A - A1; 142 Z2 = Z - Z1; 140 Z2 = Z - Z1; 143 if (A2 < 1 || Z2 < 0 || Z2 > A2) { 141 if (A2 < 1 || Z2 < 0 || Z2 > A2) { 144 FragmentsExcitationEnergy = -1.0; 142 FragmentsExcitationEnergy = -1.0; 145 continue; 143 continue; 146 } 144 } 147 M2 = G4NucleiProperties::GetNuclearMass(A2 145 M2 = G4NucleiProperties::GetNuclearMass(A2, Z2); 148 // Maximal Kinetic Energy (available energ 146 // Maximal Kinetic Energy (available energy for fragments) 149 G4double Tmax = M + U - M1 - M2 - pcorr; 147 G4double Tmax = M + U - M1 - M2 - pcorr; 150 148 151 // Check that fragment masses are less or 149 // Check that fragment masses are less or equal than total energy 152 if (Tmax < 0.0) { 150 if (Tmax < 0.0) { 153 FragmentsExcitationEnergy = -1.0; 151 FragmentsExcitationEnergy = -1.0; 154 continue; 152 continue; 155 } 153 } 156 154 157 FragmentsKineticEnergy = FissionKineticEne 155 FragmentsKineticEnergy = FissionKineticEnergy( A , Z, 158 A1, Z1, 156 A1, Z1, 159 A2, Z2, 157 A2, Z2, 160 U , Tmax); 158 U , Tmax); 161 159 162 // Excitation Energy 160 // Excitation Energy 163 // FragmentsExcitationEnergy = Tmax - Frag 161 // FragmentsExcitationEnergy = Tmax - FragmentsKineticEnergy; 164 // JMQ 04/03/09 BUG FIXED: in order to ful 162 // JMQ 04/03/09 BUG FIXED: in order to fulfill energy conservation the 165 // fragments carry the fission pairing ene 163 // fragments carry the fission pairing energy in form of 166 // excitation energy 164 // excitation energy 167 165 168 FragmentsExcitationEnergy = 166 FragmentsExcitationEnergy = 169 Tmax - FragmentsKineticEnergy + pcorr; 167 Tmax - FragmentsKineticEnergy + pcorr; 170 168 171 // Loop checking, 05-Aug-2015, Vladimir Iv 169 // Loop checking, 05-Aug-2015, Vladimir Ivanchenko 172 } while (FragmentsExcitationEnergy < 0.0 && 170 } while (FragmentsExcitationEnergy < 0.0 && ++Trials < 100); 173 171 174 if (FragmentsExcitationEnergy <= 0.0) { 172 if (FragmentsExcitationEnergy <= 0.0) { 175 throw G4HadronicException(__FILE__, __LINE 173 throw G4HadronicException(__FILE__, __LINE__, 176 "G4CompetitiveFission::BreakItUp: Excita 174 "G4CompetitiveFission::BreakItUp: Excitation energy for fragments < 0.0!"); 177 } 175 } 178 176 179 // Fragment 1 177 // Fragment 1 180 M1 += FragmentsExcitationEnergy * A1/static_ 178 M1 += FragmentsExcitationEnergy * A1/static_cast<G4double>(A); 181 // Fragment 2 179 // Fragment 2 182 M2 += FragmentsExcitationEnergy * A2/static_ 180 M2 += FragmentsExcitationEnergy * A2/static_cast<G4double>(A); 183 // primary 181 // primary 184 M += U; 182 M += U; 185 183 186 G4double etot1 = ((M - M2)*(M + M2) + M1*M1) 184 G4double etot1 = ((M - M2)*(M + M2) + M1*M1)/(2*M); 187 G4ParticleMomentum Momentum1 = 185 G4ParticleMomentum Momentum1 = 188 std::sqrt((etot1 - M1)*(etot1+M1))*G4Rando 186 std::sqrt((etot1 - M1)*(etot1+M1))*G4RandomDirection(); 189 G4LorentzVector FourMomentum1(Momentum1, eto 187 G4LorentzVector FourMomentum1(Momentum1, etot1); 190 FourMomentum1.boost(theNucleusMomentum.boost 188 FourMomentum1.boost(theNucleusMomentum.boostVector()); 191 189 192 // Create Fragments 190 // Create Fragments 193 Fragment1 = new G4Fragment( A1, Z1, FourMome 191 Fragment1 = new G4Fragment( A1, Z1, FourMomentum1); 194 if (Fragment1 != nullptr) { Fragment1->SetCr 192 if (Fragment1 != nullptr) { Fragment1->SetCreatorModelID(theSecID); } 195 theNucleusMomentum -= FourMomentum1; 193 theNucleusMomentum -= FourMomentum1; 196 theNucleus->SetZandA_asInt(Z2, A2); 194 theNucleus->SetZandA_asInt(Z2, A2); 197 theNucleus->SetMomentum(theNucleusMomentum); 195 theNucleus->SetMomentum(theNucleusMomentum); 198 theNucleus->SetCreatorModelID(theSecID); 196 theNucleus->SetCreatorModelID(theSecID); 199 return Fragment1; 197 return Fragment1; 200 } 198 } 201 199 202 G4int 200 G4int 203 G4CompetitiveFission::FissionAtomicNumber(G4in 201 G4CompetitiveFission::FissionAtomicNumber(G4int A) 204 // Calculates the atomic number of a fission 202 // Calculates the atomic number of a fission product 205 { 203 { 206 204 207 // For Simplicity reading code 205 // For Simplicity reading code 208 G4int A1 = theParam.GetA1(); 206 G4int A1 = theParam.GetA1(); 209 G4int A2 = theParam.GetA2(); 207 G4int A2 = theParam.GetA2(); 210 G4double As = theParam.GetAs(); 208 G4double As = theParam.GetAs(); 211 G4double Sigma2 = theParam.GetSigma2(); 209 G4double Sigma2 = theParam.GetSigma2(); 212 G4double SigmaS = theParam.GetSigmaS(); 210 G4double SigmaS = theParam.GetSigmaS(); 213 G4double w = theParam.GetW(); 211 G4double w = theParam.GetW(); 214 212 215 G4double C2A = A2 + 3.72*Sigma2; 213 G4double C2A = A2 + 3.72*Sigma2; 216 G4double C2S = As + 3.72*SigmaS; 214 G4double C2S = As + 3.72*SigmaS; 217 215 218 G4double C2 = 0.0; 216 G4double C2 = 0.0; 219 if (w > 1000.0 ) { C2 = C2S; } 217 if (w > 1000.0 ) { C2 = C2S; } 220 else if (w < 0.001) { C2 = C2A; } 218 else if (w < 0.001) { C2 = C2A; } 221 else { C2 = std::max(C2A,C2S 219 else { C2 = std::max(C2A,C2S); } 222 220 223 G4double C1 = A-C2; 221 G4double C1 = A-C2; 224 if (C1 < 30.0) { 222 if (C1 < 30.0) { 225 C2 = A-30.0; 223 C2 = A-30.0; 226 C1 = 30.0; 224 C1 = 30.0; 227 } 225 } 228 226 229 G4double Am1 = (As + A1)*0.5; 227 G4double Am1 = (As + A1)*0.5; 230 G4double Am2 = (A1 + A2)*0.5; 228 G4double Am2 = (A1 + A2)*0.5; 231 229 232 // Get Mass distributions as sum of symmetri 230 // Get Mass distributions as sum of symmetric and asymmetric Gasussians 233 G4double Mass1 = MassDistribution(As,A); 231 G4double Mass1 = MassDistribution(As,A); 234 G4double Mass2 = MassDistribution(Am1,A); 232 G4double Mass2 = MassDistribution(Am1,A); 235 G4double Mass3 = MassDistribution(G4double(A 233 G4double Mass3 = MassDistribution(G4double(A1),A); 236 G4double Mass4 = MassDistribution(Am2,A); 234 G4double Mass4 = MassDistribution(Am2,A); 237 G4double Mass5 = MassDistribution(G4double(A 235 G4double Mass5 = MassDistribution(G4double(A2),A); 238 // get maximal value among Mass1,...,Mass5 236 // get maximal value among Mass1,...,Mass5 239 G4double MassMax = Mass1; 237 G4double MassMax = Mass1; 240 if (Mass2 > MassMax) { MassMax = Mass2; } 238 if (Mass2 > MassMax) { MassMax = Mass2; } 241 if (Mass3 > MassMax) { MassMax = Mass3; } 239 if (Mass3 > MassMax) { MassMax = Mass3; } 242 if (Mass4 > MassMax) { MassMax = Mass4; } 240 if (Mass4 > MassMax) { MassMax = Mass4; } 243 if (Mass5 > MassMax) { MassMax = Mass5; } 241 if (Mass5 > MassMax) { MassMax = Mass5; } 244 242 245 // Sample a fragment mass number, which lies 243 // Sample a fragment mass number, which lies between C1 and C2 246 G4double xm; 244 G4double xm; 247 G4double Pm; 245 G4double Pm; 248 do { 246 do { 249 xm = C1+G4UniformRand()*(C2-C1); 247 xm = C1+G4UniformRand()*(C2-C1); 250 Pm = MassDistribution(xm,A); 248 Pm = MassDistribution(xm,A); 251 // Loop checking, 05-Aug-2015, Vladimir Iv 249 // Loop checking, 05-Aug-2015, Vladimir Ivanchenko 252 } while (MassMax*G4UniformRand() > Pm); 250 } while (MassMax*G4UniformRand() > Pm); 253 G4int ires = G4lrint(xm); 251 G4int ires = G4lrint(xm); 254 252 255 return ires; 253 return ires; 256 } 254 } 257 255 258 G4double 256 G4double 259 G4CompetitiveFission::MassDistribution(G4doubl 257 G4CompetitiveFission::MassDistribution(G4double x, G4int A) 260 // This method gives mass distribution F(x) 258 // This method gives mass distribution F(x) = F_{asym}(x)+w*F_{sym}(x) 261 // which consist of symmetric and asymmetric 259 // which consist of symmetric and asymmetric sum of gaussians components. 262 { 260 { 263 G4double y0 = (x-theParam.GetAs())/theParam. 261 G4double y0 = (x-theParam.GetAs())/theParam.GetSigmaS(); 264 G4double Xsym = LocalExp(y0); 262 G4double Xsym = LocalExp(y0); 265 263 266 G4double y1 = (x - theParam.GetA1())/thePara 264 G4double y1 = (x - theParam.GetA1())/theParam.GetSigma1(); 267 G4double y2 = (x - theParam.GetA2())/thePara 265 G4double y2 = (x - theParam.GetA2())/theParam.GetSigma2(); 268 G4double z1 = (x - A + theParam.GetA1())/the 266 G4double z1 = (x - A + theParam.GetA1())/theParam.GetSigma1(); 269 G4double z2 = (x - A + theParam.GetA2())/the 267 G4double z2 = (x - A + theParam.GetA2())/theParam.GetSigma2(); 270 G4double Xasym = LocalExp(y1) + LocalExp(y2) 268 G4double Xasym = LocalExp(y1) + LocalExp(y2) 271 + 0.5*(LocalExp(z1) + LocalExp(z2)); 269 + 0.5*(LocalExp(z1) + LocalExp(z2)); 272 270 273 G4double res; 271 G4double res; 274 G4double w = theParam.GetW(); 272 G4double w = theParam.GetW(); 275 if (w > 1000) { res = Xsym; } 273 if (w > 1000) { res = Xsym; } 276 else if (w < 0.001) { res = Xasym; } 274 else if (w < 0.001) { res = Xasym; } 277 else { res = w*Xsym+Xasym; } 275 else { res = w*Xsym+Xasym; } 278 return res; 276 return res; 279 } 277 } 280 278 281 G4int G4CompetitiveFission::FissionCharge(G4in 279 G4int G4CompetitiveFission::FissionCharge(G4int A, G4int Z, G4double Af) 282 // Calculates the charge of a fission produc 280 // Calculates the charge of a fission product for a given atomic number Af 283 { 281 { 284 static const G4double sigma = 0.6; 282 static const G4double sigma = 0.6; 285 G4double DeltaZ = 0.0; 283 G4double DeltaZ = 0.0; 286 if (Af >= 134.0) { DeltaZ = -0.45; 284 if (Af >= 134.0) { DeltaZ = -0.45; } 287 else if (Af <= (A-134.0)) { DeltaZ = 0.45; } 285 else if (Af <= (A-134.0)) { DeltaZ = 0.45; } 288 else { DeltaZ = -0.45*( 286 else { DeltaZ = -0.45*(Af-A*0.5)/(134.0-A*0.5); } 289 287 290 G4double Zmean = (Af/A)*Z + DeltaZ; 288 G4double Zmean = (Af/A)*Z + DeltaZ; 291 289 292 G4double theZ; 290 G4double theZ; 293 do { 291 do { 294 theZ = G4RandGauss::shoot(Zmean,sigma); 292 theZ = G4RandGauss::shoot(Zmean,sigma); 295 // Loop checking, 05-Aug-2015, Vladimir Iv 293 // Loop checking, 05-Aug-2015, Vladimir Ivanchenko 296 } while (theZ < 1.0 || theZ > (Z-1.0) || th 294 } while (theZ < 1.0 || theZ > (Z-1.0) || theZ > Af); 297 295 298 return G4lrint(theZ); 296 return G4lrint(theZ); 299 } 297 } 300 298 301 G4double 299 G4double 302 G4CompetitiveFission::FissionKineticEnergy(G4i 300 G4CompetitiveFission::FissionKineticEnergy(G4int A, G4int Z, 303 G4int Af1, G4int /*Zf1*/, 301 G4int Af1, G4int /*Zf1*/, 304 G4int Af2, G4int /*Zf2*/, 302 G4int Af2, G4int /*Zf2*/, 305 G4double /*U*/, G4double Tmax) 303 G4double /*U*/, G4double Tmax) 306 // Gives the kinetic energy of fission produ 304 // Gives the kinetic energy of fission products 307 { 305 { 308 // Find maximal value of A for fragments 306 // Find maximal value of A for fragments 309 G4int AfMax = std::max(Af1,Af2); 307 G4int AfMax = std::max(Af1,Af2); 310 308 311 // Weights for symmetric and asymmetric comp 309 // Weights for symmetric and asymmetric components 312 G4double Pas = 0.0; 310 G4double Pas = 0.0; 313 if (theParam.GetW() <= 1000) { 311 if (theParam.GetW() <= 1000) { 314 G4double x1 = (AfMax-theParam.GetA1())/the 312 G4double x1 = (AfMax-theParam.GetA1())/theParam.GetSigma1(); 315 G4double x2 = (AfMax-theParam.GetA2())/the 313 G4double x2 = (AfMax-theParam.GetA2())/theParam.GetSigma2(); 316 Pas = 0.5*LocalExp(x1) + LocalExp(x2); 314 Pas = 0.5*LocalExp(x1) + LocalExp(x2); 317 } 315 } 318 316 319 G4double Ps = 0.0; 317 G4double Ps = 0.0; 320 if (theParam.GetW() >= 0.001) { 318 if (theParam.GetW() >= 0.001) { 321 G4double xs = (AfMax-theParam.GetAs())/the 319 G4double xs = (AfMax-theParam.GetAs())/theParam.GetSigmaS(); 322 Ps = theParam.GetW()*LocalExp(xs); 320 Ps = theParam.GetW()*LocalExp(xs); 323 } 321 } 324 G4double Psy = (Pas + Ps > 0.0) ? Ps/(Pas+Ps 322 G4double Psy = (Pas + Ps > 0.0) ? Ps/(Pas+Ps) : 0.5; 325 323 326 // Fission fractions Xsy and Xas formed in s 324 // Fission fractions Xsy and Xas formed in symmetric and asymmetric modes 327 G4double PPas = theParam.GetSigma1() + 2.0 * 325 G4double PPas = theParam.GetSigma1() + 2.0 * theParam.GetSigma2(); 328 G4double PPsy = theParam.GetW() * theParam.G 326 G4double PPsy = theParam.GetW() * theParam.GetSigmaS(); 329 G4double Xas = (PPas + PPsy > 0.0) ? PPas/(P 327 G4double Xas = (PPas + PPsy > 0.0) ? PPas/(PPas+PPsy) : 0.5; 330 G4double Xsy = 1.0 - Xas; 328 G4double Xsy = 1.0 - Xas; 331 329 332 // Average kinetic energy for symmetric and 330 // Average kinetic energy for symmetric and asymmetric components 333 G4double Eaverage = (0.1071*(Z*Z)/G4Pow::Get 331 G4double Eaverage = (0.1071*(Z*Z)/G4Pow::GetInstance()->Z13(A) + 22.2)*CLHEP::MeV; 334 332 335 // Compute maximal average kinetic energy of 333 // Compute maximal average kinetic energy of fragments and Energy Dispersion 336 G4double TaverageAfMax; 334 G4double TaverageAfMax; 337 G4double ESigma = 10*CLHEP::MeV; 335 G4double ESigma = 10*CLHEP::MeV; 338 // Select randomly fission mode (symmetric o 336 // Select randomly fission mode (symmetric or asymmetric) 339 if (G4UniformRand() > Psy) { // Asymmetric M 337 if (G4UniformRand() > Psy) { // Asymmetric Mode 340 G4double A11 = theParam.GetA1()-0.7979*the 338 G4double A11 = theParam.GetA1()-0.7979*theParam.GetSigma1(); 341 G4double A12 = theParam.GetA1()+0.7979*the 339 G4double A12 = theParam.GetA1()+0.7979*theParam.GetSigma1(); 342 G4double A21 = theParam.GetA2()-0.7979*the 340 G4double A21 = theParam.GetA2()-0.7979*theParam.GetSigma2(); 343 G4double A22 = theParam.GetA2()+0.7979*the 341 G4double A22 = theParam.GetA2()+0.7979*theParam.GetSigma2(); 344 // scale factor 342 // scale factor 345 G4double ScaleFactor = 0.5*theParam.GetSig 343 G4double ScaleFactor = 0.5*theParam.GetSigma1()* 346 (AsymmetricRatio(A,A11)+AsymmetricRatio( 344 (AsymmetricRatio(A,A11)+AsymmetricRatio(A,A12))+ 347 theParam.GetSigma2()*(AsymmetricRatio(A, 345 theParam.GetSigma2()*(AsymmetricRatio(A,A21)+AsymmetricRatio(A,A22)); 348 // Compute average kinetic energy for frag 346 // Compute average kinetic energy for fragment with AfMax 349 TaverageAfMax = (Eaverage + 12.5 * Xsy) * 347 TaverageAfMax = (Eaverage + 12.5 * Xsy) * (PPas/ScaleFactor) * 350 AsymmetricRatio(A,G4double(AfMax)); 348 AsymmetricRatio(A,G4double(AfMax)); 351 349 352 } else { // Symmetric Mode 350 } else { // Symmetric Mode 353 G4double As0 = theParam.GetAs() + 0.7979*t 351 G4double As0 = theParam.GetAs() + 0.7979*theParam.GetSigmaS(); 354 // Compute average kinetic energy for frag 352 // Compute average kinetic energy for fragment with AfMax 355 TaverageAfMax = (Eaverage - 12.5*CLHEP::Me 353 TaverageAfMax = (Eaverage - 12.5*CLHEP::MeV*Xas) 356 *SymmetricRatio(A, G4double(AfMax))/Symm 354 *SymmetricRatio(A, G4double(AfMax))/SymmetricRatio(A, As0); 357 ESigma = 8.0*CLHEP::MeV; 355 ESigma = 8.0*CLHEP::MeV; 358 } 356 } 359 357 360 // Select randomly, in accordance with Gauss 358 // Select randomly, in accordance with Gaussian distribution, 361 // fragment kinetic energy 359 // fragment kinetic energy 362 G4double KineticEnergy; 360 G4double KineticEnergy; 363 G4int i = 0; 361 G4int i = 0; 364 do { 362 do { 365 KineticEnergy = G4RandGauss::shoot(Taverag 363 KineticEnergy = G4RandGauss::shoot(TaverageAfMax, ESigma); 366 if (++i > 100) return Eaverage; 364 if (++i > 100) return Eaverage; 367 // Loop checking, 05-Aug-2015, Vladimir Iv 365 // Loop checking, 05-Aug-2015, Vladimir Ivanchenko 368 } while (KineticEnergy < Eaverage-3.72*ESigm 366 } while (KineticEnergy < Eaverage-3.72*ESigma || 369 KineticEnergy > Eaverage+3.72*ESigma || 367 KineticEnergy > Eaverage+3.72*ESigma || 370 KineticEnergy > Tmax); 368 KineticEnergy > Tmax); 371 369 372 return KineticEnergy; 370 return KineticEnergy; 373 } 371 } 374 372 375 void G4CompetitiveFission::SetFissionBarrier(G 373 void G4CompetitiveFission::SetFissionBarrier(G4VFissionBarrier * aBarrier) 376 { 374 { 377 if (myOwnFissionBarrier) delete theFissionBa 375 if (myOwnFissionBarrier) delete theFissionBarrierPtr; 378 theFissionBarrierPtr = aBarrier; 376 theFissionBarrierPtr = aBarrier; 379 myOwnFissionBarrier = false; 377 myOwnFissionBarrier = false; 380 } 378 } 381 379 382 void 380 void 383 G4CompetitiveFission::SetEmissionStrategy(G4VE 381 G4CompetitiveFission::SetEmissionStrategy(G4VEmissionProbability * aFissionProb) 384 { 382 { 385 if (myOwnFissionProbability) delete theFissi 383 if (myOwnFissionProbability) delete theFissionProbabilityPtr; 386 theFissionProbabilityPtr = aFissionProb; 384 theFissionProbabilityPtr = aFissionProb; 387 myOwnFissionProbability = false; 385 myOwnFissionProbability = false; 388 } 386 } 389 387 390 void 388 void 391 G4CompetitiveFission::SetLevelDensityParameter 389 G4CompetitiveFission::SetLevelDensityParameter(G4VLevelDensityParameter* aLevelDensity) 392 { 390 { 393 if (myOwnLevelDensity) delete theLevelDensit 391 if (myOwnLevelDensity) delete theLevelDensityPtr; 394 theLevelDensityPtr = aLevelDensity; 392 theLevelDensityPtr = aLevelDensity; 395 myOwnLevelDensity = false; 393 myOwnLevelDensity = false; 396 } 394 } 397 395 398 396