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