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

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

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