Geant4  v4-10.4-release
 모두 클래스 네임스페이스들 파일들 함수 변수 타입정의 열거형 타입 열거형 멤버 Friends 매크로 그룹들 페이지들
G4CompetitiveFission.cc
이 파일의 문서화 페이지로 가기
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 107060 2017-11-01 15:00:04Z 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 "G4NuclearLevelData.hh"
41 #include "G4Pow.hh"
42 #include "G4Exp.hh"
43 #include "G4RandomDirection.hh"
44 #include "G4PhysicalConstants.hh"
45 
47 {
49  MyOwnFissionBarrier = true;
50 
53 
55  MyOwnLevelDensity = true;
56 
57  MaximalKineticEnergy = -1000.0;
58  FissionBarrier = 0.0;
59  FissionProbability = 0.0;
62 }
63 
65 {
69 }
70 
72 {
73  G4int anA = fragment->GetA_asInt();
74  G4int aZ = fragment->GetZ_asInt();
75  G4double ExEnergy = fragment->GetExcitationEnergy() -
77 
78  // Saddle point excitation energy ---> A = 65
79  // Fission is excluded for A < 65
80  if (anA >= 65 && ExEnergy > 0.0) {
84  theLevelDensityPtr->LevelDensityParameter(anA,aZ,ExEnergy);
88  }
89  else {
90  MaximalKineticEnergy = -1000.0;
92  FissionProbability = 0.0;
93  }
94  return FissionProbability;
95 }
96 
98 {
99  G4Fragment * Fragment1 = 0;
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();
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 = theNucleus->GetMomentum();
115 
116  // Calculate fission parameters
117  theParam.DefineParameters(A, Z, U-pcorr, FissionBarrier);
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);
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  }
148  // Maximal Kinetic Energy (available energy for fragments)
149  G4double Tmax = M + U - M1 - M2 - pcorr;
150 
151  // Check that fragment masses are less or equal than total energy
152  if (Tmax < 0.0) {
153  FragmentsExcitationEnergy = -1.0;
154  continue;
155  }
156 
157  FragmentsKineticEnergy = FissionKineticEnergy( A , Z,
158  A1, Z1,
159  A2, Z2,
160  U , Tmax);
161 
162  // Excitation Energy
163  // FragmentsExcitationEnergy = Tmax - FragmentsKineticEnergy;
164  // JMQ 04/03/09 BUG FIXED: in order to fulfill energy conservation the
165  // fragments carry the fission pairing energy in form of
166  // excitation energy
167 
168  FragmentsExcitationEnergy =
169  // Tmax - FragmentsKineticEnergy;
170  Tmax - FragmentsKineticEnergy + pcorr;
171 
172  // Loop checking, 05-Aug-2015, Vladimir Ivanchenko
173  } while (FragmentsExcitationEnergy < 0.0
174  && ++Trials < 100);
175 
176  if (FragmentsExcitationEnergy <= 0.0) {
177  throw G4HadronicException(__FILE__, __LINE__,
178  "G4CompetitiveFission::BreakItUp: Excitation energy for fragments < 0.0!");
179  }
180 
181  // Fragment 1
182  M1 += FragmentsExcitationEnergy * A1/static_cast<G4double>(A);
183  // Fragment 2
184  M2 += FragmentsExcitationEnergy * A2/static_cast<G4double>(A);
185  // primary
186  M += U;
187 
188  G4double etot1 = ((M - M2)*(M + M2) + M1*M1)/(2*M);
189  G4ParticleMomentum Momentum1 =
190  std::sqrt((etot1 - M1)*(etot1+M1))*G4RandomDirection();
191  G4LorentzVector FourMomentum1(Momentum1, etot1);
192  FourMomentum1.boost(theNucleusMomentum.boostVector());
193 
194  // Create Fragments
195  Fragment1 = new G4Fragment( A1, Z1, FourMomentum1);
196  theNucleusMomentum -= FourMomentum1;
197  theNucleus->SetZandA_asInt(Z2, A2);
198  theNucleus->SetMomentum(theNucleusMomentum);
199  return Fragment1;
200 }
201 
202 G4int
204  // Calculates the atomic number of a fission product
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 symmetric and asymmetric Gasussians
233  G4double Mass1 = MassDistribution(As,A);
234  G4double Mass2 = MassDistribution(Am1,A);
235  G4double Mass3 = MassDistribution(G4double(A1),A);
236  G4double Mass4 = MassDistribution(Am2,A);
237  G4double Mass5 = MassDistribution(G4double(A2),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 between C1 and C2
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 Ivanchenko
252  } while (MassMax*G4UniformRand() > Pm);
253  G4int ires = G4lrint(xm);
254 
255  return ires;
256 }
257 
258 G4double
260  // This method gives mass distribution F(x) = F_{asym}(x)+w*F_{sym}(x)
261  // which consist of symmetric and asymmetric sum of gaussians components.
262 {
264  G4double Xsym = G4Exp(-0.5*y0*y0);
265 
268  G4double z1 = (x - A + theParam.GetA1())/theParam.GetSigma1();
269  G4double z2 = (x - A + theParam.GetA2())/theParam.GetSigma2();
270  G4double Xasym = G4Exp(-0.5*y1*y1) + G4Exp(-0.5*y2*y2)
271  + 0.5*( G4Exp(-0.5*z1*z1) + G4Exp(-0.5*z2*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 
282  // Calculates the charge of a fission product for a given atomic number Af
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*(Af-A*0.5)/(134.0-A*0.5); }
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 Ivanchenko
296  } while (theZ < 1.0 || theZ > (Z-1.0) || theZ > Af);
297 
298  return G4lrint(theZ);
299 }
300 
301 G4double
303  G4int Af1, G4int /*Zf1*/,
304  G4int Af2, G4int /*Zf2*/,
305  G4double /*U*/, G4double Tmax)
306  // Gives the kinetic energy of fission products
307 {
308  // Find maximal value of A for fragments
309  G4int AfMax = std::max(Af1,Af2);
310 
311  // Weights for symmetric and asymmetric components
312  G4double Pas = 0.0;
313  if (theParam.GetW() <= 1000) {
314  G4double x1 = (AfMax-theParam.GetA1())/theParam.GetSigma1();
315  G4double x2 = (AfMax-theParam.GetA2())/theParam.GetSigma2();
316  Pas = 0.5*G4Exp(-0.5*x1*x1) + G4Exp(-0.5*x2*x2);
317  }
318 
319  G4double Ps = 0.0;
320  if (theParam.GetW() >= 0.001) {
321  G4double xs = (AfMax-theParam.GetAs())/theParam.GetSigmaS();
322  Ps = theParam.GetW()*G4Exp(-0.5*xs*xs);
323  }
324  G4double Psy = Ps/(Pas+Ps);
325 
326  // Fission fractions Xsy and Xas formed in symmetric and asymmetric modes
327  G4double PPas = theParam.GetSigma1() + 2.0 * theParam.GetSigma2();
328  G4double PPsy = theParam.GetW() * theParam.GetSigmaS();
329  G4double Xas = PPas / (PPas+PPsy);
330  G4double Xsy = PPsy / (PPas+PPsy);
331 
332  // Average kinetic energy for symmetric and asymmetric components
333  G4double Eaverage = (0.1071*(Z*Z)/G4Pow::GetInstance()->Z13(A) + 22.2)*CLHEP::MeV;
334 
335  // Compute maximal average kinetic energy of fragments and Energy Dispersion
336  G4double TaverageAfMax;
337  G4double ESigma = 10*CLHEP::MeV;
338  // Select randomly fission mode (symmetric or asymmetric)
339  if (G4UniformRand() > Psy) { // Asymmetric Mode
344  // scale factor
345  G4double ScaleFactor = 0.5*theParam.GetSigma1()*
346  (AsymmetricRatio(A,A11)+AsymmetricRatio(A,A12))+
348  // Compute average kinetic energy for fragment with AfMax
349  TaverageAfMax = (Eaverage + 12.5 * Xsy) * (PPas/ScaleFactor) *
350  AsymmetricRatio(A,G4double(AfMax));
351 
352  } else { // Symmetric Mode
353  G4double As0 = theParam.GetAs() + 0.7979*theParam.GetSigmaS();
354  // Compute average kinetic energy for fragment with AfMax
355  TaverageAfMax = (Eaverage - 12.5*CLHEP::MeV*Xas)
356  *SymmetricRatio(A, G4double(AfMax))/SymmetricRatio(A, As0);
357  ESigma = 8.0*CLHEP::MeV;
358  }
359 
360  // Select randomly, in accordance with Gaussian distribution,
361  // fragment kinetic energy
362  G4double KineticEnergy;
363  G4int i = 0;
364  do {
365  KineticEnergy = G4RandGauss::shoot(TaverageAfMax, ESigma);
366  if (++i > 100) return Eaverage;
367  // Loop checking, 05-Aug-2015, Vladimir Ivanchenko
368  } while (KineticEnergy < Eaverage-3.72*ESigma ||
369  KineticEnergy > Eaverage+3.72*ESigma ||
370  KineticEnergy > Tmax);
371 
372  return KineticEnergy;
373 }
374 
375 
Float_t x
Definition: compare.C:6
G4double AsymmetricRatio(G4int A, G4double A11)
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
T max(const T t1, const T t2)
brief Return the largest of the two arguments
const G4LorentzVector & GetMomentum() const
Definition: G4Fragment.hh:300
G4double GetFissionPairingCorrection(G4int A, G4int Z) const
G4double Z13(G4int Z) const
Definition: G4Pow.hh:132
G4double GetSigma1(void) const
Float_t y1[n_points_granero]
Definition: compare.C:5
G4double GetSigma2(void) const
virtual G4double LevelDensityParameter(G4int A, G4int Z, G4double U) const =0
const double C1
Float_t x1[n_points_granero]
Definition: compare.C:5
void SetZandA_asInt(G4int Znew, G4int Anew)
Definition: G4Fragment.hh:269
static G4NuclearLevelData * GetInstance()
G4int FissionCharge(G4int A, G4int Z, G4double Af)
G4int GetA_asInt() const
Definition: G4Fragment.hh:259
#define A11
G4ThreeVector G4RandomDirection()
G4double GetGroundStateMass() const
Definition: G4Fragment.hh:281
G4int GetA2(void) const
G4FissionParameters theParam
Float_t y2[n_points_geant4]
Definition: compare.C:26
Float_t Z
G4double GetSigmaS(void) const
#define A21
double G4double
Definition: G4Types.hh:76
static G4Pow * GetInstance()
Definition: G4Pow.cc:57
G4VEmissionProbability * theFissionProbabilityPtr
G4double FissionKineticEnergy(G4int A, G4int Z, G4int Af1, G4int Zf1, G4int Af2, G4int Zf2, G4double U, G4double Tmax)
#define A22
static constexpr double MeV
#define G4UniformRand()
Definition: Randomize.hh:53
double A(double temperature)
G4PairingCorrection * pairingCorrection
G4VLevelDensityParameter * theLevelDensityPtr
void SetMomentum(const G4LorentzVector &value)
Definition: G4Fragment.hh:305
static G4double GetNuclearMass(const G4double A, const G4double Z)
G4double GetAs(void) const
void DefineParameters(G4int A, G4int Z, G4double ExEnergy, G4double FissionBarrier)
G4VFissionBarrier * theFissionBarrierPtr
G4int GetA1(void) const
G4int GetZ_asInt() const
Definition: G4Fragment.hh:264
ThreeVector shoot(const G4int Ap, const G4int Af)
Double_t Z2
Double_t Z1
int G4lrint(double ad)
Definition: templates.hh:151
int G4int
Definition: G4Types.hh:78
virtual G4double FissionBarrier(G4int A, G4int Z, G4double U) const =0
virtual G4double GetEmissionProbability(G4Fragment *theNucleus)
Hep3Vector boostVector() const
virtual G4Fragment * EmittedFragment(G4Fragment *theNucleus)
G4int FissionAtomicNumber(G4int A)
#define A12
G4double MassDistribution(G4double x, G4int A)
G4PairingCorrection * GetPairingCorrection()
Float_t x2[n_points_geant4]
Definition: compare.C:26
G4double GetExcitationEnergy() const
Definition: G4Fragment.hh:276
G4double GetW(void) const
G4double SymmetricRatio(G4int A, G4double A11)
HepLorentzVector & boost(double, double, double)
virtual G4double EmissionProbability(const G4Fragment &fragment, G4double anEnergy)=0
const double C2