source: trunk/source/processes/hadronic/models/de_excitation/fission/src/G4CompetitiveFission.cc @ 1347

Last change on this file since 1347 was 1347, checked in by garnier, 13 years ago

geant4 tag 9.4
File size: 19.3 KB
Line 
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,v 1.14 2010/11/17 20:22:46 vnivanch Exp $
28// GEANT4 tag $Name: geant4-09-04-ref-00 $
29//
30// Hadronic Process: Nuclear De-excitations
31// by V. Lara (Oct 1998)
32//
33// J. M. Quesada (March 2009). Bugs fixed:
34//          - Full relativistic calculation (Lorentz boosts)
35//          - Fission pairing energy is included in fragment excitation energies
36// Now Energy and momentum are conserved in fission
37
38#include "G4CompetitiveFission.hh"
39#include "G4PairingCorrection.hh"
40#include "G4ParticleMomentum.hh"
41#include "G4Pow.hh"
42
43G4CompetitiveFission::G4CompetitiveFission() : G4VEvaporationChannel("fission")
44{
45    theFissionBarrierPtr = new G4FissionBarrier;
46    MyOwnFissionBarrier = true;
47
48    theFissionProbabilityPtr = new G4FissionProbability;
49    MyOwnFissionProbability = true;
50 
51    theLevelDensityPtr = new G4FissionLevelDensityParameter;
52    MyOwnLevelDensity = true;
53
54    MaximalKineticEnergy = -1000.0*MeV;
55    FissionBarrier = 0.0;
56    FissionProbability = 0.0;
57    LevelDensityParameter = 0.0;
58}
59
60G4CompetitiveFission::~G4CompetitiveFission()
61{
62    if (MyOwnFissionBarrier) delete theFissionBarrierPtr;
63
64    if (MyOwnFissionProbability) delete theFissionProbabilityPtr;
65
66    if (MyOwnLevelDensity) delete theLevelDensityPtr;
67}
68
69void G4CompetitiveFission::Initialize(const G4Fragment & fragment)
70{
71  G4int anA = fragment.GetA_asInt();
72  G4int aZ  = fragment.GetZ_asInt();
73  G4double ExEnergy = fragment.GetExcitationEnergy() - 
74    G4PairingCorrection::GetInstance()->GetFissionPairingCorrection(anA,aZ);
75 
76
77  // Saddle point excitation energy ---> A = 65
78  // Fission is excluded for A < 65
79  if (anA >= 65 && ExEnergy > 0.0) {
80    FissionBarrier = theFissionBarrierPtr->FissionBarrier(anA,aZ,ExEnergy);
81    MaximalKineticEnergy = ExEnergy - FissionBarrier;
82    LevelDensityParameter = 
83      theLevelDensityPtr->LevelDensityParameter(anA,aZ,ExEnergy);
84    FissionProbability = 
85      theFissionProbabilityPtr->EmissionProbability(fragment,MaximalKineticEnergy);
86    }
87  else {
88    MaximalKineticEnergy = -1000.0*MeV;
89    LevelDensityParameter = 0.0;
90    FissionProbability = 0.0;
91  }
92}
93
94G4FragmentVector * G4CompetitiveFission::BreakUp(const G4Fragment & theNucleus)
95{
96  // Nucleus data
97  // Atomic number of nucleus
98  G4int A = theNucleus.GetA_asInt();
99  // Charge of nucleus
100  G4int Z = theNucleus.GetZ_asInt();
101  //   Excitation energy (in MeV)
102  G4double U = theNucleus.GetExcitationEnergy() - 
103    G4PairingCorrection::GetInstance()->GetFissionPairingCorrection(A,Z);
104  // Check that U > 0
105  if (U <= 0.0) {
106    G4FragmentVector * theResult = new  G4FragmentVector;
107    theResult->push_back(new G4Fragment(theNucleus));
108    return theResult;
109  }
110
111  // Atomic Mass of Nucleus (in MeV)
112  G4double M = theNucleus.GetGroundStateMass();
113
114  // Nucleus Momentum
115  G4LorentzVector theNucleusMomentum = theNucleus.GetMomentum();
116
117  // Calculate fission parameters
118  G4FissionParameters theParameters(A,Z,U,FissionBarrier);
119 
120  // First fragment
121  G4int A1 = 0;
122  G4int Z1 = 0;
123  G4double M1 = 0.0;
124
125  // Second fragment
126  G4int A2 = 0;
127  G4int Z2 = 0;
128  G4double M2 = 0.0;
129
130  G4double FragmentsExcitationEnergy = 0.0;
131  G4double FragmentsKineticEnergy = 0.0;
132
133  //JMQ 04/03/09 It will be used latter to fix the bug in energy conservation
134  G4double FissionPairingEnergy=
135    G4PairingCorrection::GetInstance()->GetFissionPairingCorrection(A,Z);
136
137  G4int Trials = 0;
138  do {
139
140    // First fragment
141    A1 = FissionAtomicNumber(A,theParameters);
142    Z1 = FissionCharge(A,Z,A1);
143    M1 = G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(Z1,A1);
144
145    // Second Fragment
146    A2 = A - A1;
147    Z2 = Z - Z1;
148    if (A2 < 1 || Z2 < 0) {
149      throw G4HadronicException(__FILE__, __LINE__, 
150        "G4CompetitiveFission::BreakUp: Can't define second fragment! ");
151    }
152    M2 = G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(Z2,A2);
153
154    // Check that fragment masses are less or equal than total energy
155    if (M1 + M2 > theNucleusMomentum.e()) {
156      throw G4HadronicException(__FILE__, __LINE__, 
157        "G4CompetitiveFission::BreakUp: Fragments Mass > Total Energy");
158    }
159    // Maximal Kinetic Energy (available energy for fragments)
160    G4double Tmax = M + U - M1 - M2;
161
162    FragmentsKineticEnergy = FissionKineticEnergy( A , Z,
163                                                   A1, Z1,
164                                                   A2, Z2,
165                                                   U , Tmax,
166                                                   theParameters);
167   
168    // Excitation Energy
169    //  FragmentsExcitationEnergy = Tmax - FragmentsKineticEnergy;
170    // JMQ 04/03/09 BUG FIXED: in order to fulfill energy conservation the
171    // fragments carry the fission pairing energy in form of
172    //excitation energy
173
174    FragmentsExcitationEnergy = 
175      Tmax - FragmentsKineticEnergy+FissionPairingEnergy;
176
177  } while (FragmentsExcitationEnergy < 0.0 && Trials++ < 100);
178   
179  if (FragmentsExcitationEnergy <= 0.0) { 
180    throw G4HadronicException(__FILE__, __LINE__, 
181      "G4CompetitiveFission::BreakItUp: Excitation energy for fragments < 0.0!");
182  }
183
184  // while (FragmentsExcitationEnergy < 0 && Trials < 100);
185 
186  // Fragment 1
187  G4double U1 = FragmentsExcitationEnergy * A1/static_cast<G4double>(A);
188    // Fragment 2
189  G4double U2 = FragmentsExcitationEnergy * A2/static_cast<G4double>(A);
190
191  //JMQ  04/03/09 Full relativistic calculation is performed
192  //
193  G4double Fragment1KineticEnergy=
194    (FragmentsKineticEnergy*(FragmentsKineticEnergy+2*(M2+U2)))
195    /(2*(M1+U1+M2+U2+FragmentsKineticEnergy));
196  G4ParticleMomentum Momentum1(IsotropicVector(std::sqrt(Fragment1KineticEnergy*(Fragment1KineticEnergy+2*(M1+U1)))));
197  G4ParticleMomentum Momentum2(-Momentum1);
198  G4LorentzVector FourMomentum1(Momentum1,std::sqrt(Momentum1.mag2()+(M1+U1)*(M1+U1)));
199  G4LorentzVector FourMomentum2(Momentum2,std::sqrt(Momentum2.mag2()+(M2+U2)*(M2+U2)));
200
201  //JMQ 04/03/09 now we do Lorentz boosts (instead of Galileo boosts)
202  FourMomentum1.boost(theNucleusMomentum.boostVector());
203  FourMomentum2.boost(theNucleusMomentum.boostVector());
204   
205  //////////JMQ 04/03: Old version calculation is commented
206  // There was vioation of energy momentum conservation
207
208  //    G4double Pmax = std::sqrt( 2 * ( ( (M1+U1)*(M2+U2) ) /
209  //                            ( (M1+U1)+(M2+U2) ) ) * FragmentsKineticEnergy);
210
211  //G4ParticleMomentum momentum1 = IsotropicVector( Pmax );
212  //  G4ParticleMomentum momentum2( -momentum1 );
213
214  // Perform a Galileo boost for fragments
215  //    momentum1 += (theNucleusMomentum.boostVector() * (M1+U1));
216  //    momentum2 += (theNucleusMomentum.boostVector() * (M2+U2));
217
218
219  // Create 4-momentum for first fragment
220  // Warning!! Energy conservation is broken
221  //JMQ 04/03/09 ...NOT ANY MORE!! BUGS FIXED: Energy and momentum are NOW conserved
222  //    G4LorentzVector FourMomentum1( momentum1 , std::sqrt(momentum1.mag2() + (M1+U1)*(M1+U1)));
223
224  // Create 4-momentum for second fragment
225  // Warning!! Energy conservation is broken
226  //JMQ 04/03/09 ...NOT ANY MORE!! BUGS FIXED: Energy and momentum are NOW conserved
227  //    G4LorentzVector FourMomentum2( momentum2 , std::sqrt(momentum2.mag2() + (M2+U2)*(M2+U2)));
228
229  //////////
230
231  // Create Fragments
232  G4Fragment * Fragment1 = new G4Fragment( A1, Z1, FourMomentum1);
233  G4Fragment * Fragment2 = new G4Fragment( A2, Z2, FourMomentum2);
234
235  // Create Fragment Vector
236  G4FragmentVector * theResult = new G4FragmentVector;
237
238  theResult->push_back(Fragment1);
239  theResult->push_back(Fragment2);
240
241#ifdef debug
242  CheckConservation(theNucleus,theResult);
243#endif
244
245  return theResult;
246}
247
248G4int
249G4CompetitiveFission::FissionAtomicNumber(G4int A, 
250                                          const G4FissionParameters & theParam)
251  // Calculates the atomic number of a fission product
252{
253
254  // For Simplicity reading code
255  const G4double A1 = theParam.GetA1();
256  const G4double A2 = theParam.GetA2();
257  const G4double As = theParam.GetAs();
258  //    const G4double Sigma1 = theParam.GetSigma1();
259  const G4double Sigma2 = theParam.GetSigma2();
260  const G4double SigmaS = theParam.GetSigmaS();
261  const G4double w = theParam.GetW();
262 
263  //    G4double FasymAsym = 2.0*std::exp(-((A2-As)*(A2-As))/(2.0*Sigma2*Sigma2)) +
264  //    std::exp(-((A1-As)*(A1-As))/(2.0*Sigma1*Sigma1));
265
266  //    G4double FsymA1A2 = std::exp(-((As-(A1+A2))*(As-(A1+A2)))/(2.0*SigmaS*SigmaS));
267
268  G4double C2A = A2 + 3.72*Sigma2;
269  G4double C2S = As + 3.72*SigmaS;
270 
271  G4double C2 = 0.0;
272  if (w > 1000.0 ) C2 = C2S;
273  else if (w < 0.001) C2 = C2A;
274  else C2 =  std::max(C2A,C2S);
275
276  G4double C1 = A-C2;
277  if (C1 < 30.0) {
278    C2 = A-30.0;
279    C1 = 30.0;
280  }
281
282  G4double Am1 = (As + A1)/2.0;
283  G4double Am2 = (A1 + A2)/2.0;
284
285  // Get Mass distributions as sum of symmetric and asymmetric Gasussians
286  G4double Mass1 = MassDistribution(As,A,theParam); 
287  G4double Mass2 = MassDistribution(Am1,A,theParam); 
288  G4double Mass3 = MassDistribution(A1,A,theParam); 
289  G4double Mass4 = MassDistribution(Am2,A,theParam); 
290  G4double Mass5 = MassDistribution(A2,A,theParam); 
291  // get maximal value among Mass1,...,Mass5
292  G4double MassMax = Mass1;
293  if (Mass2 > MassMax) MassMax = Mass2;
294  if (Mass3 > MassMax) MassMax = Mass3;
295  if (Mass4 > MassMax) MassMax = Mass4;
296  if (Mass5 > MassMax) MassMax = Mass5;
297
298  // Sample a fragment mass number, which lies between C1 and C2
299  G4double m;
300  G4double Pm;
301  do {
302    m = C1+G4UniformRand()*(C2-C1);
303    Pm = MassDistribution(m,A,theParam); 
304  } while (MassMax*G4UniformRand() > Pm);
305
306  return static_cast<G4int>(m+0.5);
307}
308
309G4double
310G4CompetitiveFission::MassDistribution(G4double x, G4double A, 
311                                       const G4FissionParameters & theParam)
312  // This method gives mass distribution F(x) = F_{asym}(x)+w*F_{sym}(x)
313  // which consist of symmetric and asymmetric sum of gaussians components.
314{
315  G4double Xsym = std::exp(-0.5*(x-theParam.GetAs())*(x-theParam.GetAs())/
316                           (theParam.GetSigmaS()*theParam.GetSigmaS()));
317
318  G4double Xasym = std::exp(-0.5*(x-theParam.GetA2())*(x-theParam.GetA2())/
319                            (theParam.GetSigma2()*theParam.GetSigma2())) + 
320    std::exp(-0.5*(x-(A-theParam.GetA2()))*(x-(A-theParam.GetA2()))/
321             (theParam.GetSigma2()*theParam.GetSigma2())) +
322    0.5*std::exp(-0.5*(x-theParam.GetA1())*(x-theParam.GetA1())/
323                 (theParam.GetSigma1()*theParam.GetSigma1())) +
324    0.5*std::exp(-0.5*(x-(A-theParam.GetA1()))*(x-(A-theParam.GetA1()))/
325                 (theParam.GetSigma1()*theParam.GetSigma1()));
326
327  if (theParam.GetW() > 1000) return Xsym;
328  else if (theParam.GetW() < 0.001) return Xasym;
329  else return theParam.GetW()*Xsym+Xasym;
330}
331
332G4int G4CompetitiveFission::FissionCharge(G4double A, G4double Z,
333                                          G4double Af)
334  // Calculates the charge of a fission product for a given atomic number Af
335{
336  const G4double sigma = 0.6;
337  G4double DeltaZ = 0.0;
338  if (Af >= 134.0) DeltaZ = -0.45;                    //                      134 <= Af
339  else if (Af <= (A-134.0)) DeltaZ = 0.45;             // Af <= (A-134)
340  else DeltaZ = -0.45*(Af-(A/2.0))/(134.0-(A/2.0));   //       (A-134) < Af < 134
341
342  G4double Zmean = (Af/A)*Z + DeltaZ;
343 
344  G4double theZ;
345  do {
346    theZ = G4RandGauss::shoot(Zmean,sigma);
347  } while (theZ  < 1.0 || theZ > (Z-1.0) || theZ > Af);
348  //  return static_cast<G4int>(theZ+0.5);
349  return static_cast<G4int>(theZ+0.5);
350}
351
352G4double
353G4CompetitiveFission::FissionKineticEnergy(G4int A, G4int Z,
354                                           G4double Af1, G4double /*Zf1*/,
355                                           G4double Af2, G4double /*Zf2*/,
356                                           G4double /*U*/, G4double Tmax,
357                                           const G4FissionParameters & theParam)
358  // Gives the kinetic energy of fission products
359{
360  // Find maximal value of A for fragments
361  G4double AfMax = std::max(Af1,Af2);
362  if (AfMax < (A/2.0)) AfMax = A - AfMax;
363
364  // Weights for symmetric and asymmetric components
365  G4double Pas;
366  if (theParam.GetW() > 1000) Pas = 0.0;
367  else {
368    G4double P1 = 0.5*std::exp(-0.5*(AfMax-theParam.GetA1())*(AfMax-theParam.GetA1())/
369                               (theParam.GetSigma1()*theParam.GetSigma1()));
370
371    G4double P2 = std::exp(-0.5*(AfMax-theParam.GetA2())*(AfMax-theParam.GetA2())/
372                           (theParam.GetSigma2()*theParam.GetSigma2()));
373
374    Pas = P1+P2;
375  }
376
377  G4double Ps;
378  if (theParam.GetW() < 0.001) Ps = 0.0;
379  else {
380    Ps = theParam.GetW()*std::exp(-0.5*(AfMax-theParam.GetAs())*(AfMax-theParam.GetAs())/
381                                  (theParam.GetSigmaS()*theParam.GetSigmaS()));
382  }
383  G4double Psy = Ps/(Pas+Ps);
384
385  // Fission fractions Xsy and Xas formed in symmetric and asymmetric modes
386  G4double PPas = theParam.GetSigma1() + 2.0 * theParam.GetSigma2();
387  G4double PPsy = theParam.GetW() * theParam.GetSigmaS();
388  G4double Xas = PPas / (PPas+PPsy);
389  G4double Xsy = PPsy / (PPas+PPsy);
390
391  // Average kinetic energy for symmetric and asymmetric components
392  G4double Eaverage = 0.1071*MeV*(Z*Z)/G4Pow::GetInstance()->Z13(A) + 22.2*MeV;
393
394
395  // Compute maximal average kinetic energy of fragments and Energy Dispersion (sqrt)
396  G4double TaverageAfMax;
397  G4double ESigma;
398  // Select randomly fission mode (symmetric or asymmetric)
399  if (G4UniformRand() > Psy) { // Asymmetric Mode
400    G4double A11 = theParam.GetA1()-0.7979*theParam.GetSigma1();
401    G4double A12 = theParam.GetA1()+0.7979*theParam.GetSigma1();
402    G4double A21 = theParam.GetA2()-0.7979*theParam.GetSigma2();
403    G4double A22 = theParam.GetA2()+0.7979*theParam.GetSigma2();
404    // scale factor
405    G4double ScaleFactor = 0.5*theParam.GetSigma1()*(AsymmetricRatio(A,A11)+AsymmetricRatio(A,A12))+
406      theParam.GetSigma2()*(AsymmetricRatio(A,A21)+AsymmetricRatio(A,A22));
407    // Compute average kinetic energy for fragment with AfMax
408    TaverageAfMax = (Eaverage + 12.5 * Xsy) * (PPas/ScaleFactor) * AsymmetricRatio(A,AfMax);
409    ESigma = 10.0*MeV; // MeV
410
411  } else { // Symmetric Mode
412    G4double As0 = theParam.GetAs() + 0.7979*theParam.GetSigmaS();
413    // scale factor
414    G4double ScaleFactor = theParam.GetW()*theParam.GetSigmaS()*SymmetricRatio(A,As0);
415    // Compute average kinetic energy for fragment with AfMax
416    TaverageAfMax = (Eaverage - 12.5*MeV*Xas) * (PPsy/ScaleFactor) * SymmetricRatio(A,AfMax);
417    ESigma = 8.0*MeV;
418  }
419
420
421  // Select randomly, in accordance with Gaussian distribution, fragment kinetic energy
422  G4double KineticEnergy;
423  G4int i = 0;
424  do {
425    KineticEnergy = G4RandGauss::shoot(TaverageAfMax,ESigma);
426    if (i++ > 100) return Eaverage;
427  } while (KineticEnergy < Eaverage-3.72*ESigma || 
428           KineticEnergy > Eaverage+3.72*ESigma ||
429           KineticEnergy > Tmax);
430 
431  return KineticEnergy;
432}
433
434G4double G4CompetitiveFission::AsymmetricRatio(G4int A, G4double A11)
435{
436  const G4double B1 = 23.5;
437  const G4double A00 = 134.0;
438  return Ratio(G4double(A),A11,B1,A00);
439}
440
441G4double G4CompetitiveFission::SymmetricRatio(G4int A, G4double A11)
442{
443  const G4double B1 = 5.32;
444  const G4double A00 = A/2.0;
445  return Ratio(G4double(A),A11,B1,A00);
446}
447
448G4double G4CompetitiveFission::Ratio(G4double A, G4double A11,
449                                     G4double B1, G4double A00) 
450{
451  if (A == 0.0) {
452    throw G4HadronicException(__FILE__, __LINE__, 
453                              "G4CompetitiveFission::Ratio: A == 0!");
454  }
455  if (A11 >= A/2.0 && A11 <= (A00+10.0)) {
456    return 1.0-B1*((A11-A00)/A)*((A11-A00)/A);
457  } else {
458    return 1.0-B1*(10.0/A)*(10.0/A)-2.0*(10.0/A)*B1*((A11-A00-10.0)/A);
459  }
460}
461
462G4ThreeVector G4CompetitiveFission::IsotropicVector(const G4double Magnitude)
463  // Samples a isotropic random vectorwith a magnitud given by Magnitude.
464  // By default Magnitude = 1.0
465{
466  G4double CosTheta = 1.0 - 2.0*G4UniformRand();
467  G4double SinTheta = std::sqrt(1.0 - CosTheta*CosTheta);
468  G4double Phi = twopi*G4UniformRand();
469  G4ThreeVector Vector(Magnitude*std::cos(Phi)*SinTheta,
470                       Magnitude*std::sin(Phi)*SinTheta,
471                       Magnitude*CosTheta);
472  return Vector;
473}
474
475#ifdef debug
476void G4CompetitiveFission::CheckConservation(const G4Fragment & theInitialState,
477                                             G4FragmentVector * Result) const
478{
479    G4double ProductsEnergy =0;
480    G4ThreeVector ProductsMomentum;
481    G4int ProductsA = 0;
482    G4int ProductsZ = 0;
483    G4FragmentVector::iterator h;
484    for (h = Result->begin(); h != Result->end(); h++) {
485        G4LorentzVector tmp = (*h)->GetMomentum();
486        ProductsEnergy += tmp.e();
487        ProductsMomentum += tmp.vect();
488        ProductsA += static_cast<G4int>((*h)->GetA());
489        ProductsZ += static_cast<G4int>((*h)->GetZ());
490    }
491
492    if (ProductsA != theInitialState.GetA()) {
493        G4cout << "!!!!!!!!!! Baryonic Number Conservation Violation !!!!!!!!!!" << G4endl;
494        G4cout << "G4CompetitiveFission.cc: Barionic Number Conservation test for fission fragments" 
495               << G4endl; 
496        G4cout << "Initial A = " << theInitialState.GetA() 
497               << "   Fragments A = " << ProductsA << "   Diference --> " 
498               << theInitialState.GetA() - ProductsA << G4endl;
499    }
500    if (ProductsZ != theInitialState.GetZ()) {
501        G4cout << "!!!!!!!!!! Charge Conservation Violation !!!!!!!!!!" << G4endl;
502        G4cout << "G4CompetitiveFission.cc: Charge Conservation test for fission fragments" 
503               << G4endl; 
504        G4cout << "Initial Z = " << theInitialState.GetZ() 
505               << "   Fragments Z = " << ProductsZ << "   Diference --> " 
506               << theInitialState.GetZ() - ProductsZ << G4endl;
507    }
508    if (std::abs(ProductsEnergy-theInitialState.GetMomentum().e()) > 1.0*keV) {
509        G4cout << "!!!!!!!!!! Energy Conservation Violation !!!!!!!!!!" << G4endl;
510        G4cout << "G4CompetitiveFission.cc: Energy Conservation test for fission fragments" 
511               << G4endl; 
512        G4cout << "Initial E = " << theInitialState.GetMomentum().e()/MeV << " MeV"
513               << "   Fragments E = " << ProductsEnergy/MeV  << " MeV   Diference --> " 
514               << (theInitialState.GetMomentum().e() - ProductsEnergy)/MeV << " MeV" << G4endl;
515    } 
516    if (std::abs(ProductsMomentum.x()-theInitialState.GetMomentum().x()) > 1.0*keV || 
517        std::abs(ProductsMomentum.y()-theInitialState.GetMomentum().y()) > 1.0*keV ||
518        std::abs(ProductsMomentum.z()-theInitialState.GetMomentum().z()) > 1.0*keV) {
519        G4cout << "!!!!!!!!!! Momentum Conservation Violation !!!!!!!!!!" << G4endl;
520        G4cout << "G4CompetitiveFission.cc: Momentum Conservation test for fission fragments" 
521               << G4endl; 
522        G4cout << "Initial P = " << theInitialState.GetMomentum().vect() << " MeV"
523               << "   Fragments P = " << ProductsMomentum  << " MeV   Diference --> " 
524               << theInitialState.GetMomentum().vect() - ProductsMomentum << " MeV" << G4endl;
525    }
526    return;
527}
528#endif
529
530
531
532
Note: See TracBrowser for help on using the repository browser.