Geant4  10.01.p03
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 85841 2014-11-05 15:35:06Z 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 "G4PhysicalConstants.hh"
42 #include "G4SystemOfUnits.hh"
43 
45 {
47  MyOwnFissionBarrier = true;
48 
51 
53  MyOwnLevelDensity = true;
54 
55  MaximalKineticEnergy = -1000.0*MeV;
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 
77  // Saddle point excitation energy ---> A = 65
78  // Fission is excluded for A < 65
79  if (anA >= 65 && ExEnergy > 0.0) {
83  theLevelDensityPtr->LevelDensityParameter(anA,aZ,ExEnergy);
86  }
87  else {
88  MaximalKineticEnergy = -1000.0*MeV;
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  G4FissionParameters theParameters(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  //JMQ 04/03/09 It will be used latter to fix the bug in energy conservation
141  G4double FissionPairingEnergy=
143 
144  G4int Trials = 0;
145  do {
146 
147  // First fragment
148  A1 = FissionAtomicNumber(A,theParameters);
149  Z1 = FissionCharge(A,Z,A1);
151 
152  // Second Fragment
153  A2 = A - A1;
154  Z2 = Z - Z1;
155  if (A2 < 1 || Z2 < 0 || Z2 > A2) {
156  FragmentsExcitationEnergy = -1.0;
157  continue;
158  }
160  // Maximal Kinetic Energy (available energy for fragments)
161  G4double Tmax = M + U - M1 - M2;
162 
163  // Check that fragment masses are less or equal than total energy
164  if (Tmax < 0.0) {
165  FragmentsExcitationEnergy = -1.0;
166  continue;
167  }
168 
169  FragmentsKineticEnergy = FissionKineticEnergy( A , Z,
170  A1, Z1,
171  A2, Z2,
172  U , Tmax,
173  theParameters);
174 
175  // Excitation Energy
176  // FragmentsExcitationEnergy = Tmax - FragmentsKineticEnergy;
177  // JMQ 04/03/09 BUG FIXED: in order to fulfill energy conservation the
178  // fragments carry the fission pairing energy in form of
179  // excitation energy
180 
181  FragmentsExcitationEnergy =
182  Tmax - FragmentsKineticEnergy+FissionPairingEnergy;
183 
184  } while (FragmentsExcitationEnergy < 0.0 && Trials++ < 100);
185 
186  if (FragmentsExcitationEnergy <= 0.0) {
187  throw G4HadronicException(__FILE__, __LINE__,
188  "G4CompetitiveFission::BreakItUp: Excitation energy for fragments < 0.0!");
189  }
190 
191  // Fragment 1
192  M1 += FragmentsExcitationEnergy * A1/static_cast<G4double>(A);
193  // Fragment 2
194  M2 += FragmentsExcitationEnergy * A2/static_cast<G4double>(A);
195  // primary
196  M += U;
197 
198  G4double etot1 = (M*M - M2*M2 + M1*M1)/(2*M);
199  G4ParticleMomentum Momentum1(IsotropicVector(std::sqrt((etot1 - M1)*(etot1+M1))));
200  G4LorentzVector FourMomentum1(Momentum1, etot1);
201  FourMomentum1.boost(theNucleusMomentum.boostVector());
202 
203  // Create Fragments
204  Fragment1 = new G4Fragment( A1, Z1, FourMomentum1);
205  theNucleusMomentum -= FourMomentum1;
206  theNucleus->SetZandA_asInt(Z2, A2);
207  theNucleus->SetMomentum(theNucleusMomentum);
208  return Fragment1;
209 }
210 
211 G4int
213  const G4FissionParameters & theParam)
214  // Calculates the atomic number of a fission product
215 {
216 
217  // For Simplicity reading code
218  G4double A1 = theParam.GetA1();
219  G4double A2 = theParam.GetA2();
220  G4double As = theParam.GetAs();
221  G4double Sigma2 = theParam.GetSigma2();
222  G4double SigmaS = theParam.GetSigmaS();
223  G4double w = theParam.GetW();
224 
225  // G4double FasymAsym = 2.0*std::exp(-((A2-As)*(A2-As))/(2.0*Sigma2*Sigma2)) +
226  // std::exp(-((A1-As)*(A1-As))/(2.0*Sigma1*Sigma1));
227 
228  // G4double FsymA1A2 = std::exp(-((As-(A1+A2))*(As-(A1+A2)))/(2.0*SigmaS*SigmaS));
229 
230  G4double C2A = A2 + 3.72*Sigma2;
231  G4double C2S = As + 3.72*SigmaS;
232 
233  G4double C2 = 0.0;
234  if (w > 1000.0 ) C2 = C2S;
235  else if (w < 0.001) C2 = C2A;
236  else C2 = std::max(C2A,C2S);
237 
238  G4double C1 = A-C2;
239  if (C1 < 30.0) {
240  C2 = A-30.0;
241  C1 = 30.0;
242  }
243 
244  G4double Am1 = (As + A1)/2.0;
245  G4double Am2 = (A1 + A2)/2.0;
246 
247  // Get Mass distributions as sum of symmetric and asymmetric Gasussians
248  G4double Mass1 = MassDistribution(As,A,theParam);
249  G4double Mass2 = MassDistribution(Am1,A,theParam);
250  G4double Mass3 = MassDistribution(A1,A,theParam);
251  G4double Mass4 = MassDistribution(Am2,A,theParam);
252  G4double Mass5 = MassDistribution(A2,A,theParam);
253  // get maximal value among Mass1,...,Mass5
254  G4double MassMax = Mass1;
255  if (Mass2 > MassMax) MassMax = Mass2;
256  if (Mass3 > MassMax) MassMax = Mass3;
257  if (Mass4 > MassMax) MassMax = Mass4;
258  if (Mass5 > MassMax) MassMax = Mass5;
259 
260  // Sample a fragment mass number, which lies between C1 and C2
261  G4double xm;
262  G4double Pm;
263  do {
264  xm = C1+G4UniformRand()*(C2-C1);
265  Pm = MassDistribution(xm,A,theParam);
266  } while (MassMax*G4UniformRand() > Pm);
267  G4int ires = G4lrint(xm);
268 
269  return ires;
270 }
271 
272 G4double
274  const G4FissionParameters & theParam)
275  // This method gives mass distribution F(x) = F_{asym}(x)+w*F_{sym}(x)
276  // which consist of symmetric and asymmetric sum of gaussians components.
277 {
278  G4double Xsym = std::exp(-0.5*(x-theParam.GetAs())*(x-theParam.GetAs())/
279  (theParam.GetSigmaS()*theParam.GetSigmaS()));
280 
281  G4double Xasym = std::exp(-0.5*(x-theParam.GetA2())*(x-theParam.GetA2())/
282  (theParam.GetSigma2()*theParam.GetSigma2())) +
283  std::exp(-0.5*(x-(A-theParam.GetA2()))*(x-(A-theParam.GetA2()))/
284  (theParam.GetSigma2()*theParam.GetSigma2())) +
285  0.5*std::exp(-0.5*(x-theParam.GetA1())*(x-theParam.GetA1())/
286  (theParam.GetSigma1()*theParam.GetSigma1())) +
287  0.5*std::exp(-0.5*(x-(A-theParam.GetA1()))*(x-(A-theParam.GetA1()))/
288  (theParam.GetSigma1()*theParam.GetSigma1()));
289 
290  if (theParam.GetW() > 1000) return Xsym;
291  else if (theParam.GetW() < 0.001) return Xasym;
292  else return theParam.GetW()*Xsym+Xasym;
293 }
294 
296  G4double Af)
297  // Calculates the charge of a fission product for a given atomic number Af
298 {
299  static const G4double sigma = 0.6;
300  G4double DeltaZ = 0.0;
301  if (Af >= 134.0) DeltaZ = -0.45; // 134 <= Af
302  else if (Af <= (A-134.0)) DeltaZ = 0.45; // Af <= (A-134)
303  else DeltaZ = -0.45*(Af-(A/2.0))/(134.0-(A/2.0)); // (A-134) < Af < 134
304 
305  G4double Zmean = (Af/A)*Z + DeltaZ;
306 
307  G4double theZ;
308  do {
309  theZ = G4RandGauss::shoot(Zmean,sigma);
310  } while (theZ < 1.0 || theZ > (Z-1.0) || theZ > Af);
311  // return static_cast<G4int>(theZ+0.5);
312  return static_cast<G4int>(theZ+0.5);
313 }
314 
315 G4double
317  G4double Af1, G4double /*Zf1*/,
318  G4double Af2, G4double /*Zf2*/,
319  G4double /*U*/, G4double Tmax,
320  const G4FissionParameters & theParam)
321  // Gives the kinetic energy of fission products
322 {
323  // Find maximal value of A for fragments
324  G4double AfMax = std::max(Af1,Af2);
325  if (AfMax < (A/2.0)) AfMax = A - AfMax;
326 
327  // Weights for symmetric and asymmetric components
328  G4double Pas;
329  if (theParam.GetW() > 1000) Pas = 0.0;
330  else {
331  G4double P1 = 0.5*std::exp(-0.5*(AfMax-theParam.GetA1())*(AfMax-theParam.GetA1())/
332  (theParam.GetSigma1()*theParam.GetSigma1()));
333 
334  G4double P2 = std::exp(-0.5*(AfMax-theParam.GetA2())*(AfMax-theParam.GetA2())/
335  (theParam.GetSigma2()*theParam.GetSigma2()));
336 
337  Pas = P1+P2;
338  }
339 
340  G4double Ps;
341  if (theParam.GetW() < 0.001) Ps = 0.0;
342  else {
343  Ps = theParam.GetW()*std::exp(-0.5*(AfMax-theParam.GetAs())*(AfMax-theParam.GetAs())/
344  (theParam.GetSigmaS()*theParam.GetSigmaS()));
345  }
346  G4double Psy = Ps/(Pas+Ps);
347 
348  // Fission fractions Xsy and Xas formed in symmetric and asymmetric modes
349  G4double PPas = theParam.GetSigma1() + 2.0 * theParam.GetSigma2();
350  G4double PPsy = theParam.GetW() * theParam.GetSigmaS();
351  G4double Xas = PPas / (PPas+PPsy);
352  G4double Xsy = PPsy / (PPas+PPsy);
353 
354  // Average kinetic energy for symmetric and asymmetric components
355  G4double Eaverage = 0.1071*MeV*(Z*Z)/G4Pow::GetInstance()->Z13(A) + 22.2*MeV;
356 
357  // Compute maximal average kinetic energy of fragments and Energy Dispersion (sqrt)
358  G4double TaverageAfMax;
359  G4double ESigma = 10*MeV;
360  // Select randomly fission mode (symmetric or asymmetric)
361  if (G4UniformRand() > Psy) { // Asymmetric Mode
362  G4double A11 = theParam.GetA1()-0.7979*theParam.GetSigma1();
363  G4double A12 = theParam.GetA1()+0.7979*theParam.GetSigma1();
364  G4double A21 = theParam.GetA2()-0.7979*theParam.GetSigma2();
365  G4double A22 = theParam.GetA2()+0.7979*theParam.GetSigma2();
366  // scale factor
367  G4double ScaleFactor = 0.5*theParam.GetSigma1()*
368  (AsymmetricRatio(A,A11)+AsymmetricRatio(A,A12))+
369  theParam.GetSigma2()*(AsymmetricRatio(A,A21)+AsymmetricRatio(A,A22));
370  // Compute average kinetic energy for fragment with AfMax
371  TaverageAfMax = (Eaverage + 12.5 * Xsy) * (PPas/ScaleFactor) *
372  AsymmetricRatio(A,AfMax);
373 
374  } else { // Symmetric Mode
375  G4double As0 = theParam.GetAs() + 0.7979*theParam.GetSigmaS();
376  // scale factor
377  G4double ScaleFactor = theParam.GetW()*theParam.GetSigmaS()*SymmetricRatio(A,As0);
378  // Compute average kinetic energy for fragment with AfMax
379  TaverageAfMax = (Eaverage - 12.5*MeV*Xas) * (PPsy/ScaleFactor) *
380  SymmetricRatio(A,AfMax);
381  ESigma = 8.0*MeV;
382  }
383 
384 
385  // Select randomly, in accordance with Gaussian distribution, fragment kinetic energy
386  G4double KineticEnergy;
387  G4int i = 0;
388  do {
389  KineticEnergy = G4RandGauss::shoot(TaverageAfMax,ESigma);
390  if (i++ > 100) return Eaverage;
391  } while (KineticEnergy < Eaverage-3.72*ESigma ||
392  KineticEnergy > Eaverage+3.72*ESigma ||
393  KineticEnergy > Tmax);
394 
395  return KineticEnergy;
396 }
397 
399 {
400  static const G4double B1 = 23.5;
401  static const G4double A00 = 134.0;
402  return Ratio(G4double(A),A11,B1,A00);
403 }
404 
406 {
407  static const G4double B1 = 5.32;
408  const G4double A00 = A/2.0;
409  return Ratio(G4double(A),A11,B1,A00);
410 }
411 
413  G4double B1, G4double A00)
414 {
415  if (A == 0.0) {
416  throw G4HadronicException(__FILE__, __LINE__,
417  "G4CompetitiveFission::Ratio: A == 0!");
418  }
419  if (A11 >= A/2.0 && A11 <= (A00+10.0)) {
420  return 1.0-B1*((A11-A00)/A)*((A11-A00)/A);
421  } else {
422  return 1.0-B1*(10.0/A)*(10.0/A)-2.0*(10.0/A)*B1*((A11-A00-10.0)/A);
423  }
424 }
425 
427  // Samples a isotropic random vectorwith a magnitud given by Magnitude.
428  // By default Magnitude = 1.0
429 {
430  G4double CosTheta = 1.0 - 2.0*G4UniformRand();
431  G4double SinTheta = std::sqrt(1.0 - CosTheta*CosTheta);
432  G4double Phi = twopi*G4UniformRand();
433  G4ThreeVector Vector(Magnitude*std::cos(Phi)*SinTheta,
434  Magnitude*std::sin(Phi)*SinTheta,
435  Magnitude*CosTheta);
436  return Vector;
437 }
438 
439 #ifdef debug
440 void G4CompetitiveFission::CheckConservation(const G4Fragment & theInitialState,
441  G4FragmentVector * Result) const
442 {
443  G4double ProductsEnergy =0;
444  G4ThreeVector ProductsMomentum;
445  G4int ProductsA = 0;
446  G4int ProductsZ = 0;
447  G4FragmentVector::iterator h;
448  for (h = Result->begin(); h != Result->end(); h++) {
449  G4LorentzVector tmp = (*h)->GetMomentum();
450  ProductsEnergy += tmp.e();
451  ProductsMomentum += tmp.vect();
452  ProductsA += (*h)->GetA_asInt();
453  ProductsZ += (*h)->GetZ_asInt();
454  }
455 
456  if (ProductsA != theInitialState.GetA_asInt()) {
457  G4cout << "!!!!!!!!!! Baryonic Number Conservation Violation !!!!!!!!!!"
458  << G4endl;
459  G4cout << "G4CompetitiveFission: Baryon Number Conservation test for fission fragments"
460  << G4endl;
461  G4cout << "Initial A = " << theInitialState.GetA_asInt()
462  << " Fragments A = " << ProductsA << " Diference --> "
463  << theInitialState.GetA_asInt() - ProductsA << G4endl;
464  }
465  if (ProductsZ != theInitialState.GetZ_asInt()) {
466  G4cout << "!!!!!!!!!! Charge Conservation Violation !!!!!!!!!!" << G4endl;
467  G4cout << "G4CompetitiveFission.cc: Charge Conservation test for fission fragments"
468  << G4endl;
469  G4cout << "Initial Z = " << theInitialState.GetZ_asInt()
470  << " Fragments Z = " << ProductsZ << " Diference --> "
471  << theInitialState.GetZ() - ProductsZ << G4endl;
472  }
473  if (std::fabs(ProductsEnergy-theInitialState.GetMomentum().e()) > 1.0*keV) {
474  G4cout << "!!!!!!!!!! Energy Conservation Violation !!!!!!!!!!" << G4endl;
475  G4cout << "G4CompetitiveFission.cc: Energy Conservation test for fission fragments"
476  << G4endl;
477  G4cout << "Initial E = " << theInitialState.GetMomentum().e()/MeV << " MeV"
478  << " Fragments E = " << ProductsEnergy/MeV << " MeV Diference --> "
479  << (theInitialState.GetMomentum().e() - ProductsEnergy)/MeV << " MeV" << G4endl;
480  }
481  if (std::fabs(ProductsMomentum.x()-theInitialState.GetMomentum().x()) > 1.0*keV ||
482  std::fabs(ProductsMomentum.y()-theInitialState.GetMomentum().y()) > 1.0*keV ||
483  std::fabs(ProductsMomentum.z()-theInitialState.GetMomentum().z()) > 1.0*keV) {
484  G4cout << "!!!!!!!!!! Momentum Conservation Violation !!!!!!!!!!" << G4endl;
485  G4cout << "G4CompetitiveFission.cc: Momentum Conservation test for fission fragments"
486  << G4endl;
487  G4cout << "Initial P = " << theInitialState.GetMomentum().vect() << " MeV"
488  << " Fragments P = " << ProductsMomentum << " MeV Diference --> "
489  << theInitialState.GetMomentum().vect() - ProductsMomentum << " MeV" << G4endl;
490  }
491  return;
492 }
493 #endif
494 
495 
496 
497 
#define A21
static G4Pow * GetInstance()
Definition: G4Pow.cc:55
ThreeVector shoot(const G4int Ap, const G4int Af)
G4double GetSigma2(void) const
static const double MeV
Definition: G4SIunits.hh:193
G4double AsymmetricRatio(G4int A, G4double A11)
static G4double GetNuclearMass(const G4double A, const G4double Z)
static const G4double * P1[nN]
G4double GetSigma1(void) const
CLHEP::Hep3Vector G4ThreeVector
#define A00
G4double MassDistribution(G4double x, G4double A, const G4FissionParameters &theParam)
#define A22
virtual G4double GetEmissionProbability(G4Fragment *theNucleus)
G4int FissionCharge(G4double A, G4double Z, G4double Af)
G4VEmissionProbability * theFissionProbabilityPtr
virtual G4double FissionBarrier(G4int A, G4int Z, const G4double U)=0
G4int FissionAtomicNumber(G4int A, const G4FissionParameters &theParam)
int G4int
Definition: G4Types.hh:78
virtual G4FragmentVector * BreakUp(const G4Fragment &theNucleus)
#define A12
G4PairingCorrection * pairingCorrection
#define G4UniformRand()
Definition: Randomize.hh:93
G4GLOB_DLL std::ostream G4cout
G4double Z13(G4int Z) const
Definition: G4Pow.hh:127
G4int GetA_asInt() const
Definition: G4Fragment.hh:243
G4double GetSigmaS(void) const
const G4LorentzVector & GetMomentum() const
Definition: G4Fragment.hh:276
void SetMomentum(const G4LorentzVector &value)
Definition: G4Fragment.hh:281
virtual G4Fragment * EmittedFragment(G4Fragment *theNucleus)
std::vector< G4Fragment * > G4FragmentVector
Definition: G4Fragment.hh:65
G4double Ratio(G4double A, G4double A11, G4double B1, G4double A00)
G4double SymmetricRatio(G4int A, G4double A11)
G4double GetGroundStateMass() const
Definition: G4Fragment.hh:265
G4VLevelDensityParameter * theLevelDensityPtr
static const G4double A[nN]
G4double GetA1(void) const
G4VFissionBarrier * theFissionBarrierPtr
#define C1
G4ThreeVector IsotropicVector(G4double Magnitude=1.0)
static G4PairingCorrection * GetInstance()
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:253
G4double FissionKineticEnergy(G4int A, G4int Z, G4double Af1, G4double Zf1, G4double Af2, G4double Zf2, G4double U, G4double Tmax, const G4FissionParameters &theParam)
virtual G4double LevelDensityParameter(G4int A, G4int Z, G4double U) const =0
#define A11
G4int GetZ_asInt() const
Definition: G4Fragment.hh:248
static const G4double * P2[nN]
#define G4endl
Definition: G4ios.hh:61
static const double keV
Definition: G4SIunits.hh:195
double G4double
Definition: G4Types.hh:76
G4double GetW(void) const
virtual G4double EmissionProbability(const G4Fragment &fragment, const G4double anEnergy)=0
G4ThreeVector G4ParticleMomentum
G4double GetExcitationEnergy() const
Definition: G4Fragment.hh:260
G4double GetZ() const
Definition: G4Fragment.hh:303
CLHEP::HepLorentzVector G4LorentzVector
G4double GetA2(void) const