| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 //
27 // Hadronic Process: Nuclear De-excitations
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
35 #include "G4CompetitiveFission.hh"
36 #include "G4PairingCorrection.hh"
37 #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"
45 #include "Randomize.hh"
46 #include "G4RandomDirection.hh"
47 #include "G4PhysicalConstants.hh"
48 #include "G4PhysicsModelCatalog.hh"
49
50 G4CompetitiveFission::G4CompetitiveFission() :
51 {
52 theFissionBarrierPtr = new G4FissionBarrier;
53 theFissionProbabilityPtr = new G4FissionProb
54 theLevelDensityPtr = new G4FissionLevelDensi
55 pairingCorrection = G4NuclearLevelData::GetI
56 theSecID = G4PhysicsModelCatalog::GetModelID
57 }
58
59 G4CompetitiveFission::~G4CompetitiveFission()
60 {
61 if (myOwnFissionBarrier) delete theFissionBa
62 if (myOwnFissionProbability) delete theFissi
63 if (myOwnLevelDensity) delete theLevelDensit
64 }
65
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
76 {
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 }
96
97 G4Fragment* G4CompetitiveFission::EmittedFragm
98 {
99 G4Fragment * Fragment1 = nullptr;
100 // Nucleus data
101 // Atomic number of nucleus
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
110 // Atomic Mass of Nucleus (in MeV)
111 G4double M = theNucleus->GetGroundStateMass(
112
113 // Nucleus Momentum
114 G4LorentzVector theNucleusMomentum = theNucl
115
116 // Calculate fission parameters
117 theParam.DefineParameters(A, Z, U-pcorr, fis
118
119 // First fragment
120 G4int A1 = 0;
121 G4int Z1 = 0;
122 G4double M1 = 0.0;
123
124 // Second fragment
125 G4int A2 = 0;
126 G4int Z2 = 0;
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 }
147 M2 = G4NucleiProperties::GetNuclearMass(A2
148 // Maximal Kinetic Energy (available energ
149 G4double Tmax = M + U - M1 - M2 - pcorr;
150
151 // Check that fragment masses are less or
152 if (Tmax < 0.0) {
153 FragmentsExcitationEnergy = -1.0;
154 continue;
155 }
156
157 FragmentsKineticEnergy = FissionKineticEne
158 A1, Z1,
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
168 FragmentsExcitationEnergy =
169 Tmax - FragmentsKineticEnergy + pcorr;
170
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
179 // Fragment 1
180 M1 += FragmentsExcitationEnergy * A1/static_
181 // Fragment 2
182 M2 += FragmentsExcitationEnergy * A2/static_
183 // primary
184 M += U;
185
186 G4double etot1 = ((M - M2)*(M + M2) + M1*M1)
187 G4ParticleMomentum Momentum1 =
188 std::sqrt((etot1 - M1)*(etot1+M1))*G4Rando
189 G4LorentzVector FourMomentum1(Momentum1, eto
190 FourMomentum1.boost(theNucleusMomentum.boost
191
192 // Create Fragments
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
215 G4double C2A = A2 + 3.72*Sigma2;
216 G4double C2S = As + 3.72*SigmaS;
217
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
290 G4double Zmean = (Af/A)*Z + DeltaZ;
291
292 G4double theZ;
293 do {
294 theZ = G4RandGauss::shoot(Zmean,sigma);
295 // Loop checking, 05-Aug-2015, Vladimir Iv
296 } while (theZ < 1.0 || theZ > (Z-1.0) || th
297
298 return G4lrint(theZ);
299 }
300
301 G4double
302 G4CompetitiveFission::FissionKineticEnergy(G4i
303 G4int Af1, G4int /*Zf1*/,
304 G4int Af2, G4int /*Zf2*/,
305 G4double /*U*/, G4double Tmax)
306 // Gives the kinetic energy of fission produ
307 {
308 // Find maximal value of A for fragments
309 G4int AfMax = std::max(Af1,Af2);
310
311 // Weights for symmetric and asymmetric comp
312 G4double Pas = 0.0;
313 if (theParam.GetW() <= 1000) {
314 G4double x1 = (AfMax-theParam.GetA1())/the
315 G4double x2 = (AfMax-theParam.GetA2())/the
316 Pas = 0.5*LocalExp(x1) + LocalExp(x2);
317 }
318
319 G4double Ps = 0.0;
320 if (theParam.GetW() >= 0.001) {
321 G4double xs = (AfMax-theParam.GetAs())/the
322 Ps = theParam.GetW()*LocalExp(xs);
323 }
324 G4double Psy = (Pas + Ps > 0.0) ? Ps/(Pas+Ps
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
372 return KineticEnergy;
373 }
374
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
398