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

Last change on this file since 1228 was 1228, checked in by garnier, 14 years ago

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