Geant4  10.02.p02
G4CompetitiveFission.cc
Go to the documentation of this file.
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 //
27 // $Id: G4CompetitiveFission.cc 91834 2015-08-07 07:24:22Z gcosmo $
28 //
29 // Hadronic Process: Nuclear De-excitations
30 // by V. Lara (Oct 1998)
31 //
32 // J. M. Quesada (March 2009). Bugs fixed:
33 // - Full relativistic calculation (Lorentz boosts)
34 // - Fission pairing energy is included in fragment excitation energies
35 // Now Energy and momentum are conserved in fission
36 
37 #include "G4CompetitiveFission.hh"
38 #include "G4PairingCorrection.hh"
39 #include "G4ParticleMomentum.hh"
40 #include "G4Pow.hh"
41 #include "G4Exp.hh"
42 #include "G4PhysicalConstants.hh"
43 
45 {
47  MyOwnFissionBarrier = true;
48 
51 
53  MyOwnLevelDensity = true;
54 
55  MaximalKineticEnergy = -1000.0;
56  FissionBarrier = 0.0;
57  FissionProbability = 0.0;
60 }
61 
63 {
67 }
68 
70 {
71  G4int anA = fragment->GetA_asInt();
72  G4int aZ = fragment->GetZ_asInt();
73  G4double ExEnergy = fragment->GetExcitationEnergy() -
75 
76  // Saddle point excitation energy ---> A = 65
77  // Fission is excluded for A < 65
78  if (anA >= 65 && ExEnergy > 0.0) {
82  theLevelDensityPtr->LevelDensityParameter(anA,aZ,ExEnergy);
86  }
87  else {
88  MaximalKineticEnergy = -1000.0;
90  FissionProbability = 0.0;
91  }
92  return FissionProbability;
93 }
94 
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;
103 }
104 
106 {
107  G4Fragment * Fragment1 = 0;
108  // Nucleus data
109  // Atomic number of nucleus
110  G4int A = theNucleus->GetA_asInt();
111  // Charge of nucleus
112  G4int Z = theNucleus->GetZ_asInt();
113  // Excitation energy (in MeV)
114  G4double U = theNucleus->GetExcitationEnergy();
116  if (U <= pcorr) { return Fragment1; }
117 
118  // Atomic Mass of Nucleus (in MeV)
119  G4double M = theNucleus->GetGroundStateMass();
120 
121  // Nucleus Momentum
122  G4LorentzVector theNucleusMomentum = theNucleus->GetMomentum();
123 
124  // Calculate fission parameters
125  theParam.DefineParameters(A, Z, U-pcorr, FissionBarrier);
126 
127  // First fragment
128  G4int A1 = 0;
129  G4int Z1 = 0;
130  G4double M1 = 0.0;
131 
132  // Second fragment
133  G4int A2 = 0;
134  G4int Z2 = 0;
135  G4double M2 = 0.0;
136 
137  G4double FragmentsExcitationEnergy = 0.0;
138  G4double FragmentsKineticEnergy = 0.0;
139 
140  G4int Trials = 0;
141  do {
142 
143  // First fragment
144  A1 = FissionAtomicNumber(A);
145  Z1 = FissionCharge(A, Z, A1);
147 
148  // Second Fragment
149  A2 = A - A1;
150  Z2 = Z - Z1;
151  if (A2 < 1 || Z2 < 0 || Z2 > A2) {
152  FragmentsExcitationEnergy = -1.0;
153  continue;
154  }
156  // Maximal Kinetic Energy (available energy for fragments)
157  G4double Tmax = M + U - M1 - M2 - pcorr;
158 
159  // Check that fragment masses are less or equal than total energy
160  if (Tmax < 0.0) {
161  FragmentsExcitationEnergy = -1.0;
162  continue;
163  }
164 
165  FragmentsKineticEnergy = FissionKineticEnergy( A , Z,
166  A1, Z1,
167  A2, Z2,
168  U , Tmax);
169 
170  // Excitation Energy
171  // FragmentsExcitationEnergy = Tmax - FragmentsKineticEnergy;
172  // JMQ 04/03/09 BUG FIXED: in order to fulfill energy conservation the
173  // fragments carry the fission pairing energy in form of
174  // excitation energy
175 
176  FragmentsExcitationEnergy =
177  // Tmax - FragmentsKineticEnergy;
178  Tmax - FragmentsKineticEnergy + pcorr;
179 
180  // Loop checking, 05-Aug-2015, Vladimir Ivanchenko
181  } while (FragmentsExcitationEnergy < 0.0
182  && ++Trials < 100);
183 
184  if (FragmentsExcitationEnergy <= 0.0) {
185  throw G4HadronicException(__FILE__, __LINE__,
186  "G4CompetitiveFission::BreakItUp: Excitation energy for fragments < 0.0!");
187  }
188 
189  // Fragment 1
190  M1 += FragmentsExcitationEnergy * A1/static_cast<G4double>(A);
191  // Fragment 2
192  M2 += FragmentsExcitationEnergy * A2/static_cast<G4double>(A);
193  // primary
194  M += U;
195 
196  G4double etot1 = ((M - M2)*(M + M2) + M1*M1)/(2*M);
197  G4ParticleMomentum Momentum1(IsotropicVector(std::sqrt((etot1 - M1)*(etot1+M1))));
198  G4LorentzVector FourMomentum1(Momentum1, etot1);
199  FourMomentum1.boost(theNucleusMomentum.boostVector());
200 
201  // Create Fragments
202  Fragment1 = new G4Fragment( A1, Z1, FourMomentum1);
203  theNucleusMomentum -= FourMomentum1;
204  theNucleus->SetZandA_asInt(Z2, A2);
205  theNucleus->SetMomentum(theNucleusMomentum);
206  return Fragment1;
207 }
208 
209 G4int
211  // Calculates the atomic number of a fission product
212 {
213 
214  // For Simplicity reading code
215  G4int A1 = theParam.GetA1();
216  G4int A2 = theParam.GetA2();
217  G4double As = theParam.GetAs();
218  G4double Sigma2 = theParam.GetSigma2();
219  G4double SigmaS = theParam.GetSigmaS();
220  G4double w = theParam.GetW();
221 
222  G4double C2A = A2 + 3.72*Sigma2;
223  G4double C2S = As + 3.72*SigmaS;
224 
225  G4double C2 = 0.0;
226  if (w > 1000.0 ) { C2 = C2S; }
227  else if (w < 0.001) { C2 = C2A; }
228  else { C2 = std::max(C2A,C2S); }
229 
230  G4double C1 = A-C2;
231  if (C1 < 30.0) {
232  C2 = A-30.0;
233  C1 = 30.0;
234  }
235 
236  G4double Am1 = (As + A1)*0.5;
237  G4double Am2 = (A1 + A2)*0.5;
238 
239  // Get Mass distributions as sum of symmetric and asymmetric Gasussians
240  G4double Mass1 = MassDistribution(As,A);
241  G4double Mass2 = MassDistribution(Am1,A);
242  G4double Mass3 = MassDistribution(G4double(A1),A);
243  G4double Mass4 = MassDistribution(Am2,A);
244  G4double Mass5 = MassDistribution(G4double(A2),A);
245  // get maximal value among Mass1,...,Mass5
246  G4double MassMax = Mass1;
247  if (Mass2 > MassMax) { MassMax = Mass2; }
248  if (Mass3 > MassMax) { MassMax = Mass3; }
249  if (Mass4 > MassMax) { MassMax = Mass4; }
250  if (Mass5 > MassMax) { MassMax = Mass5; }
251 
252  // Sample a fragment mass number, which lies between C1 and C2
253  G4double xm;
254  G4double Pm;
255  do {
256  xm = C1+G4UniformRand()*(C2-C1);
257  Pm = MassDistribution(xm,A);
258  // Loop checking, 05-Aug-2015, Vladimir Ivanchenko
259  } while (MassMax*G4UniformRand() > Pm);
260  G4int ires = G4lrint(xm);
261 
262  return ires;
263 }
264 
265 G4double
267  // This method gives mass distribution F(x) = F_{asym}(x)+w*F_{sym}(x)
268  // which consist of symmetric and asymmetric sum of gaussians components.
269 {
271  G4double Xsym = G4Exp(-0.5*y0*y0);
272 
273  G4double y1 = (x - theParam.GetA1())/theParam.GetSigma1();
274  G4double y2 = (x - theParam.GetA2())/theParam.GetSigma2();
275  G4double z1 = (x - A + theParam.GetA1())/theParam.GetSigma1();
276  G4double z2 = (x - A + theParam.GetA2())/theParam.GetSigma2();
277  G4double Xasym = G4Exp(-0.5*y1*y1) + G4Exp(-0.5*y2*y2)
278  + 0.5*( G4Exp(-0.5*z1*z1) + G4Exp(-0.5*z2*z2));
279 
280  G4double res;
281  G4double w = theParam.GetW();
282  if (w > 1000) { res = Xsym; }
283  else if (w < 0.001) { res = Xasym; }
284  else { res = w*Xsym+Xasym; }
285  return res;
286 }
287 
289  // Calculates the charge of a fission product for a given atomic number Af
290 {
291  static const G4double sigma = 0.6;
292  G4double DeltaZ = 0.0;
293  if (Af >= 134.0) { DeltaZ = -0.45; }
294  else if (Af <= (A-134.0)) { DeltaZ = 0.45; }
295  else { DeltaZ = -0.45*(Af-A*0.5)/(134.0-A*0.5); }
296 
297  G4double Zmean = (Af/A)*Z + DeltaZ;
298 
299  G4double theZ;
300  do {
301  theZ = G4RandGauss::shoot(Zmean,sigma);
302  // Loop checking, 05-Aug-2015, Vladimir Ivanchenko
303  } while (theZ < 1.0 || theZ > (Z-1.0) || theZ > Af);
304 
305  return G4lrint(theZ);
306 }
307 
308 G4double
310  G4int Af1, G4int /*Zf1*/,
311  G4int Af2, G4int /*Zf2*/,
312  G4double /*U*/, G4double Tmax)
313  // Gives the kinetic energy of fission products
314 {
315  // Find maximal value of A for fragments
316  G4int AfMax = std::max(Af1,Af2);
317 
318  // Weights for symmetric and asymmetric components
319  G4double Pas = 0.0;
320  if (theParam.GetW() <= 1000) {
321  G4double x1 = (AfMax-theParam.GetA1())/theParam.GetSigma1();
322  G4double x2 = (AfMax-theParam.GetA2())/theParam.GetSigma2();
323  Pas = 0.5*G4Exp(-0.5*x1*x1) + G4Exp(-0.5*x2*x2);
324  }
325 
326  G4double Ps = 0.0;
327  if (theParam.GetW() >= 0.001) {
328  G4double xs = (AfMax-theParam.GetAs())/theParam.GetSigmaS();
329  Ps = theParam.GetW()*G4Exp(-0.5*xs*xs);
330  }
331  G4double Psy = Ps/(Pas+Ps);
332 
333  // Fission fractions Xsy and Xas formed in symmetric and asymmetric modes
334  G4double PPas = theParam.GetSigma1() + 2.0 * theParam.GetSigma2();
335  G4double PPsy = theParam.GetW() * theParam.GetSigmaS();
336  G4double Xas = PPas / (PPas+PPsy);
337  G4double Xsy = PPsy / (PPas+PPsy);
338 
339  // Average kinetic energy for symmetric and asymmetric components
340  G4double Eaverage = (0.1071*(Z*Z)/G4Pow::GetInstance()->Z13(A) + 22.2)*CLHEP::MeV;
341 
342  // Compute maximal average kinetic energy of fragments and Energy Dispersion
343  G4double TaverageAfMax;
344  G4double ESigma = 10*CLHEP::MeV;
345  // Select randomly fission mode (symmetric or asymmetric)
346  if (G4UniformRand() > Psy) { // Asymmetric Mode
351  // scale factor
352  G4double ScaleFactor = 0.5*theParam.GetSigma1()*
353  (AsymmetricRatio(A,A11)+AsymmetricRatio(A,A12))+
355  // Compute average kinetic energy for fragment with AfMax
356  TaverageAfMax = (Eaverage + 12.5 * Xsy) * (PPas/ScaleFactor) *
357  AsymmetricRatio(A,G4double(AfMax));
358 
359  } else { // Symmetric Mode
360  G4double As0 = theParam.GetAs() + 0.7979*theParam.GetSigmaS();
361  // Compute average kinetic energy for fragment with AfMax
362  TaverageAfMax = (Eaverage - 12.5*CLHEP::MeV*Xas)
363  *SymmetricRatio(A, G4double(AfMax))/SymmetricRatio(A, As0);
364  ESigma = 8.0*CLHEP::MeV;
365  }
366 
367  // Select randomly, in accordance with Gaussian distribution,
368  // fragment kinetic energy
369  G4double KineticEnergy;
370  G4int i = 0;
371  do {
372  KineticEnergy = G4RandGauss::shoot(TaverageAfMax, ESigma);
373  if (++i > 100) return Eaverage;
374  // Loop checking, 05-Aug-2015, Vladimir Ivanchenko
375  } while (KineticEnergy < Eaverage-3.72*ESigma ||
376  KineticEnergy > Eaverage+3.72*ESigma ||
377  KineticEnergy > Tmax);
378 
379  return KineticEnergy;
380 }
381 
382 
#define A21
const double C2
static G4Pow * GetInstance()
Definition: G4Pow.cc:55
G4ThreeVector IsotropicVector(G4double Magnitude)
ThreeVector shoot(const G4int Ap, const G4int Af)
G4double GetSigma2(void) const
static const double MeV
Definition: G4SIunits.hh:211
G4int FissionAtomicNumber(G4int A)
G4double AsymmetricRatio(G4int A, G4double A11)
static G4double GetNuclearMass(const G4double A, const G4double Z)
const double C1
G4double GetSigma1(void) const
#define A22
virtual G4double GetEmissionProbability(G4Fragment *theNucleus)
G4VEmissionProbability * theFissionProbabilityPtr
const G4double w[NPOINTSGL]
G4FissionParameters theParam
virtual G4double FissionBarrier(G4int A, G4int Z, const G4double U)=0
int G4int
Definition: G4Types.hh:78
virtual G4FragmentVector * BreakUp(const G4Fragment &theNucleus)
#define A12
G4double MassDistribution(G4double x, G4int A)
G4PairingCorrection * pairingCorrection
G4int GetA2(void) const
#define G4UniformRand()
Definition: Randomize.hh:97
G4double Z13(G4int Z) const
Definition: G4Pow.hh:127
double A(double temperature)
G4int GetA_asInt() const
Definition: G4Fragment.hh:256
G4double GetSigmaS(void) const
const G4LorentzVector & GetMomentum() const
Definition: G4Fragment.hh:289
void SetMomentum(const G4LorentzVector &value)
Definition: G4Fragment.hh:294
virtual G4Fragment * EmittedFragment(G4Fragment *theNucleus)
std::vector< G4Fragment * > G4FragmentVector
Definition: G4Fragment.hh:63
G4double SymmetricRatio(G4int A, G4double A11)
G4double GetGroundStateMass() const
Definition: G4Fragment.hh:278
G4VLevelDensityParameter * theLevelDensityPtr
G4int FissionCharge(G4int A, G4int Z, G4double Af)
G4VFissionBarrier * theFissionBarrierPtr
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
static G4PairingCorrection * GetInstance()
G4int GetA1(void) const
G4double GetFissionPairingCorrection(G4int A, G4int Z) const
G4double GetAs(void) const
int G4lrint(double ad)
Definition: templates.hh:163
T max(const T t1, const T t2)
brief Return the largest of the two arguments
void SetZandA_asInt(G4int Znew, G4int Anew)
Definition: G4Fragment.hh:266
const G4double x[NPOINTSGL]
virtual G4double LevelDensityParameter(G4int A, G4int Z, G4double U) const =0
#define A11
G4int GetZ_asInt() const
Definition: G4Fragment.hh:261
double G4double
Definition: G4Types.hh:76
G4double GetW(void) const
virtual G4double EmissionProbability(const G4Fragment &fragment, const G4double anEnergy)=0
G4double FissionKineticEnergy(G4int A, G4int Z, G4int Af1, G4int Zf1, G4int Af2, G4int Zf2, G4double U, G4double Tmax)
G4ThreeVector G4ParticleMomentum
G4double GetExcitationEnergy() const
Definition: G4Fragment.hh:273
CLHEP::HepLorentzVector G4LorentzVector
void DefineParameters(G4int A, G4int Z, G4double ExEnergy, G4double FissionBarrier)