G4CompetitiveFission.cc

Go to the documentation of this file.
00001 //
00002 // ********************************************************************
00003 // * License and Disclaimer                                           *
00004 // *                                                                  *
00005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
00006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
00007 // * conditions of the Geant4 Software License,  included in the file *
00008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
00009 // * include a list of copyright holders.                             *
00010 // *                                                                  *
00011 // * Neither the authors of this software system, nor their employing *
00012 // * institutes,nor the agencies providing financial support for this *
00013 // * work  make  any representation or  warranty, express or implied, *
00014 // * regarding  this  software system or assume any liability for its *
00015 // * use.  Please see the license in the file  LICENSE  and URL above *
00016 // * for the full disclaimer and the limitation of liability.         *
00017 // *                                                                  *
00018 // * This  code  implementation is the result of  the  scientific and *
00019 // * technical work of the GEANT4 collaboration.                      *
00020 // * By using,  copying,  modifying or  distributing the software (or *
00021 // * any work based  on the software)  you  agree  to acknowledge its *
00022 // * use  in  resulting  scientific  publications,  and indicate your *
00023 // * acceptance of all terms of the Geant4 Software license.          *
00024 // ********************************************************************
00025 //
00026 //
00027 // $Id$
00028 //
00029 // Hadronic Process: Nuclear De-excitations
00030 // by V. Lara (Oct 1998)
00031 //
00032 // J. M. Quesada (March 2009). Bugs fixed:
00033 //          - Full relativistic calculation (Lorentz boosts)
00034 //          - Fission pairing energy is included in fragment excitation energies
00035 // Now Energy and momentum are conserved in fission 
00036 
00037 #include "G4CompetitiveFission.hh"
00038 #include "G4PairingCorrection.hh"
00039 #include "G4ParticleMomentum.hh"
00040 #include "G4Pow.hh"
00041 #include "G4PhysicalConstants.hh"
00042 #include "G4SystemOfUnits.hh"
00043 
00044 G4CompetitiveFission::G4CompetitiveFission() : G4VEvaporationChannel("fission")
00045 {
00046     theFissionBarrierPtr = new G4FissionBarrier;
00047     MyOwnFissionBarrier = true;
00048 
00049     theFissionProbabilityPtr = new G4FissionProbability;
00050     MyOwnFissionProbability = true;
00051   
00052     theLevelDensityPtr = new G4FissionLevelDensityParameter;
00053     MyOwnLevelDensity = true;
00054 
00055     MaximalKineticEnergy = -1000.0*MeV;
00056     FissionBarrier = 0.0;
00057     FissionProbability = 0.0;
00058     LevelDensityParameter = 0.0;
00059 }
00060 
00061 G4CompetitiveFission::~G4CompetitiveFission()
00062 {
00063     if (MyOwnFissionBarrier) delete theFissionBarrierPtr;
00064 
00065     if (MyOwnFissionProbability) delete theFissionProbabilityPtr;
00066 
00067     if (MyOwnLevelDensity) delete theLevelDensityPtr;
00068 }
00069 
00070 G4double G4CompetitiveFission::GetEmissionProbability(G4Fragment* fragment)
00071 {
00072   G4int anA = fragment->GetA_asInt();
00073   G4int aZ  = fragment->GetZ_asInt();
00074   G4double ExEnergy = fragment->GetExcitationEnergy() - 
00075     G4PairingCorrection::GetInstance()->GetFissionPairingCorrection(anA,aZ);
00076   
00077 
00078   // Saddle point excitation energy ---> A = 65
00079   // Fission is excluded for A < 65
00080   if (anA >= 65 && ExEnergy > 0.0) {
00081     FissionBarrier = theFissionBarrierPtr->FissionBarrier(anA,aZ,ExEnergy);
00082     MaximalKineticEnergy = ExEnergy - FissionBarrier;
00083     LevelDensityParameter = 
00084       theLevelDensityPtr->LevelDensityParameter(anA,aZ,ExEnergy);
00085     FissionProbability = 
00086       theFissionProbabilityPtr->EmissionProbability(*fragment,MaximalKineticEnergy);
00087     }
00088   else {
00089     MaximalKineticEnergy = -1000.0*MeV;
00090     LevelDensityParameter = 0.0;
00091     FissionProbability = 0.0;
00092   }
00093   return FissionProbability;
00094 }
00095 
00096 G4FragmentVector * G4CompetitiveFission::BreakUp(const G4Fragment & theNucleus)
00097 {
00098   // Nucleus data
00099   // Atomic number of nucleus
00100   G4int A = theNucleus.GetA_asInt();
00101   // Charge of nucleus
00102   G4int Z = theNucleus.GetZ_asInt();
00103   //   Excitation energy (in MeV)
00104   G4double U = theNucleus.GetExcitationEnergy() - 
00105     G4PairingCorrection::GetInstance()->GetFissionPairingCorrection(A,Z);
00106   // Check that U > 0
00107   if (U <= 0.0) {
00108     G4FragmentVector * theResult = new  G4FragmentVector;
00109     theResult->push_back(new G4Fragment(theNucleus));
00110     return theResult;
00111   }
00112 
00113   // Atomic Mass of Nucleus (in MeV)
00114   G4double M = theNucleus.GetGroundStateMass();
00115 
00116   // Nucleus Momentum
00117   G4LorentzVector theNucleusMomentum = theNucleus.GetMomentum();
00118 
00119   // Calculate fission parameters
00120   G4FissionParameters theParameters(A,Z,U,FissionBarrier);
00121   
00122   // First fragment
00123   G4int A1 = 0;
00124   G4int Z1 = 0;
00125   G4double M1 = 0.0;
00126 
00127   // Second fragment
00128   G4int A2 = 0;
00129   G4int Z2 = 0;
00130   G4double M2 = 0.0;
00131 
00132   G4double FragmentsExcitationEnergy = 0.0;
00133   G4double FragmentsKineticEnergy = 0.0;
00134 
00135   //JMQ 04/03/09 It will be used latter to fix the bug in energy conservation
00136   G4double FissionPairingEnergy=
00137     G4PairingCorrection::GetInstance()->GetFissionPairingCorrection(A,Z);
00138 
00139   G4int Trials = 0;
00140   do {
00141 
00142     // First fragment 
00143     A1 = FissionAtomicNumber(A,theParameters);
00144     Z1 = FissionCharge(A,Z,A1);
00145     M1 = G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(Z1,A1);
00146 
00147     // Second Fragment
00148     A2 = A - A1;
00149     Z2 = Z - Z1;
00150     if (A2 < 1 || Z2 < 0) {
00151       throw G4HadronicException(__FILE__, __LINE__, 
00152         "G4CompetitiveFission::BreakUp: Can't define second fragment! ");
00153     }
00154     M2 = G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(Z2,A2);
00155 
00156     // Check that fragment masses are less or equal than total energy
00157     if (M1 + M2 > theNucleusMomentum.e()) {
00158       throw G4HadronicException(__FILE__, __LINE__, 
00159         "G4CompetitiveFission::BreakUp: Fragments Mass > Total Energy");
00160     }
00161     // Maximal Kinetic Energy (available energy for fragments)
00162     G4double Tmax = M + U - M1 - M2;
00163 
00164     FragmentsKineticEnergy = FissionKineticEnergy( A , Z,
00165                                                    A1, Z1,
00166                                                    A2, Z2,
00167                                                    U , Tmax,
00168                                                    theParameters);
00169     
00170     // Excitation Energy
00171     //  FragmentsExcitationEnergy = Tmax - FragmentsKineticEnergy;
00172     // JMQ 04/03/09 BUG FIXED: in order to fulfill energy conservation the
00173     // fragments carry the fission pairing energy in form of 
00174     //excitation energy
00175 
00176     FragmentsExcitationEnergy = 
00177       Tmax - FragmentsKineticEnergy+FissionPairingEnergy;
00178 
00179   } while (FragmentsExcitationEnergy < 0.0 && Trials++ < 100);
00180     
00181   if (FragmentsExcitationEnergy <= 0.0) { 
00182     throw G4HadronicException(__FILE__, __LINE__, 
00183       "G4CompetitiveFission::BreakItUp: Excitation energy for fragments < 0.0!");
00184   }
00185 
00186   // while (FragmentsExcitationEnergy < 0 && Trials < 100);
00187   
00188   // Fragment 1
00189   G4double U1 = FragmentsExcitationEnergy * A1/static_cast<G4double>(A);
00190     // Fragment 2
00191   G4double U2 = FragmentsExcitationEnergy * A2/static_cast<G4double>(A);
00192 
00193   //JMQ  04/03/09 Full relativistic calculation is performed
00194   //
00195   G4double Fragment1KineticEnergy=
00196     (FragmentsKineticEnergy*(FragmentsKineticEnergy+2*(M2+U2)))
00197     /(2*(M1+U1+M2+U2+FragmentsKineticEnergy));
00198   G4ParticleMomentum Momentum1(IsotropicVector(std::sqrt(Fragment1KineticEnergy*(Fragment1KineticEnergy+2*(M1+U1)))));
00199   G4ParticleMomentum Momentum2(-Momentum1);
00200   G4LorentzVector FourMomentum1(Momentum1,std::sqrt(Momentum1.mag2()+(M1+U1)*(M1+U1)));
00201   G4LorentzVector FourMomentum2(Momentum2,std::sqrt(Momentum2.mag2()+(M2+U2)*(M2+U2)));
00202 
00203   //JMQ 04/03/09 now we do Lorentz boosts (instead of Galileo boosts)
00204   FourMomentum1.boost(theNucleusMomentum.boostVector());
00205   FourMomentum2.boost(theNucleusMomentum.boostVector());
00206     
00208   // There was vioation of energy momentum conservation
00209 
00210   //    G4double Pmax = std::sqrt( 2 * ( ( (M1+U1)*(M2+U2) ) /
00211   //                            ( (M1+U1)+(M2+U2) ) ) * FragmentsKineticEnergy);
00212 
00213   //G4ParticleMomentum momentum1 = IsotropicVector( Pmax );
00214   //  G4ParticleMomentum momentum2( -momentum1 );
00215 
00216   // Perform a Galileo boost for fragments
00217   //    momentum1 += (theNucleusMomentum.boostVector() * (M1+U1));
00218   //    momentum2 += (theNucleusMomentum.boostVector() * (M2+U2));
00219 
00220 
00221   // Create 4-momentum for first fragment
00222   // Warning!! Energy conservation is broken
00223   //JMQ 04/03/09 ...NOT ANY MORE!! BUGS FIXED: Energy and momentum are NOW conserved 
00224   //    G4LorentzVector FourMomentum1( momentum1 , std::sqrt(momentum1.mag2() + (M1+U1)*(M1+U1)));
00225 
00226   // Create 4-momentum for second fragment
00227   // Warning!! Energy conservation is broken
00228   //JMQ 04/03/09 ...NOT ANY MORE!! BUGS FIXED: Energy and momentum are NOW conserved
00229   //    G4LorentzVector FourMomentum2( momentum2 , std::sqrt(momentum2.mag2() + (M2+U2)*(M2+U2)));
00230 
00232 
00233   // Create Fragments
00234   G4Fragment * Fragment1 = new G4Fragment( A1, Z1, FourMomentum1);
00235   G4Fragment * Fragment2 = new G4Fragment( A2, Z2, FourMomentum2);
00236 
00237   // Create Fragment Vector
00238   G4FragmentVector * theResult = new G4FragmentVector;
00239 
00240   theResult->push_back(Fragment1);
00241   theResult->push_back(Fragment2);
00242 
00243 #ifdef debug
00244   CheckConservation(theNucleus,theResult);
00245 #endif
00246 
00247   return theResult;
00248 }
00249 
00250 G4int 
00251 G4CompetitiveFission::FissionAtomicNumber(G4int A, 
00252                                           const G4FissionParameters & theParam)
00253   // Calculates the atomic number of a fission product
00254 {
00255 
00256   // For Simplicity reading code
00257   const G4double A1 = theParam.GetA1();
00258   const G4double A2 = theParam.GetA2();
00259   const G4double As = theParam.GetAs();
00260   //    const G4double Sigma1 = theParam.GetSigma1();
00261   const G4double Sigma2 = theParam.GetSigma2();
00262   const G4double SigmaS = theParam.GetSigmaS();
00263   const G4double w = theParam.GetW();
00264   
00265   //    G4double FasymAsym = 2.0*std::exp(-((A2-As)*(A2-As))/(2.0*Sigma2*Sigma2)) + 
00266   //    std::exp(-((A1-As)*(A1-As))/(2.0*Sigma1*Sigma1));
00267 
00268   //    G4double FsymA1A2 = std::exp(-((As-(A1+A2))*(As-(A1+A2)))/(2.0*SigmaS*SigmaS));
00269 
00270   G4double C2A = A2 + 3.72*Sigma2;
00271   G4double C2S = As + 3.72*SigmaS;
00272   
00273   G4double C2 = 0.0;
00274   if (w > 1000.0 ) C2 = C2S;
00275   else if (w < 0.001) C2 = C2A;
00276   else C2 =  std::max(C2A,C2S);
00277 
00278   G4double C1 = A-C2;
00279   if (C1 < 30.0) {
00280     C2 = A-30.0;
00281     C1 = 30.0;
00282   }
00283 
00284   G4double Am1 = (As + A1)/2.0;
00285   G4double Am2 = (A1 + A2)/2.0;
00286 
00287   // Get Mass distributions as sum of symmetric and asymmetric Gasussians
00288   G4double Mass1 = MassDistribution(As,A,theParam); 
00289   G4double Mass2 = MassDistribution(Am1,A,theParam); 
00290   G4double Mass3 = MassDistribution(A1,A,theParam); 
00291   G4double Mass4 = MassDistribution(Am2,A,theParam); 
00292   G4double Mass5 = MassDistribution(A2,A,theParam); 
00293   // get maximal value among Mass1,...,Mass5
00294   G4double MassMax = Mass1;
00295   if (Mass2 > MassMax) MassMax = Mass2;
00296   if (Mass3 > MassMax) MassMax = Mass3;
00297   if (Mass4 > MassMax) MassMax = Mass4;
00298   if (Mass5 > MassMax) MassMax = Mass5;
00299 
00300   // Sample a fragment mass number, which lies between C1 and C2
00301   G4double xm;
00302   G4double Pm;
00303   do {
00304     xm = C1+G4UniformRand()*(C2-C1);
00305     Pm = MassDistribution(xm,A,theParam); 
00306   } while (MassMax*G4UniformRand() > Pm);
00307   G4int ires = G4lrint(xm);
00308 
00309   return ires;
00310 }
00311 
00312 G4double 
00313 G4CompetitiveFission::MassDistribution(G4double x, G4double A, 
00314                                        const G4FissionParameters & theParam)
00315   // This method gives mass distribution F(x) = F_{asym}(x)+w*F_{sym}(x)
00316   // which consist of symmetric and asymmetric sum of gaussians components.
00317 {
00318   G4double Xsym = std::exp(-0.5*(x-theParam.GetAs())*(x-theParam.GetAs())/
00319                            (theParam.GetSigmaS()*theParam.GetSigmaS()));
00320 
00321   G4double Xasym = std::exp(-0.5*(x-theParam.GetA2())*(x-theParam.GetA2())/
00322                             (theParam.GetSigma2()*theParam.GetSigma2())) + 
00323     std::exp(-0.5*(x-(A-theParam.GetA2()))*(x-(A-theParam.GetA2()))/
00324              (theParam.GetSigma2()*theParam.GetSigma2())) +
00325     0.5*std::exp(-0.5*(x-theParam.GetA1())*(x-theParam.GetA1())/
00326                  (theParam.GetSigma1()*theParam.GetSigma1())) +
00327     0.5*std::exp(-0.5*(x-(A-theParam.GetA1()))*(x-(A-theParam.GetA1()))/
00328                  (theParam.GetSigma1()*theParam.GetSigma1()));
00329 
00330   if (theParam.GetW() > 1000) return Xsym;
00331   else if (theParam.GetW() < 0.001) return Xasym;
00332   else return theParam.GetW()*Xsym+Xasym;
00333 }
00334 
00335 G4int G4CompetitiveFission::FissionCharge(G4double A, G4double Z,
00336                                           G4double Af)
00337   // Calculates the charge of a fission product for a given atomic number Af
00338 {
00339   const G4double sigma = 0.6;
00340   G4double DeltaZ = 0.0;
00341   if (Af >= 134.0) DeltaZ = -0.45;                    //                      134 <= Af
00342   else if (Af <= (A-134.0)) DeltaZ = 0.45;             // Af <= (A-134) 
00343   else DeltaZ = -0.45*(Af-(A/2.0))/(134.0-(A/2.0));   //       (A-134) < Af < 134
00344 
00345   G4double Zmean = (Af/A)*Z + DeltaZ;
00346  
00347   G4double theZ;
00348   do {
00349     theZ = G4RandGauss::shoot(Zmean,sigma);
00350   } while (theZ  < 1.0 || theZ > (Z-1.0) || theZ > Af);
00351   //  return static_cast<G4int>(theZ+0.5);
00352   return static_cast<G4int>(theZ+0.5);
00353 }
00354 
00355 G4double 
00356 G4CompetitiveFission::FissionKineticEnergy(G4int A, G4int Z,
00357                                            G4double Af1, G4double /*Zf1*/,
00358                                            G4double Af2, G4double /*Zf2*/,
00359                                            G4double /*U*/, G4double Tmax,
00360                                            const G4FissionParameters & theParam)
00361   // Gives the kinetic energy of fission products
00362 {
00363   // Find maximal value of A for fragments
00364   G4double AfMax = std::max(Af1,Af2);
00365   if (AfMax < (A/2.0)) AfMax = A - AfMax;
00366 
00367   // Weights for symmetric and asymmetric components
00368   G4double Pas;
00369   if (theParam.GetW() > 1000) Pas = 0.0;
00370   else {
00371     G4double P1 = 0.5*std::exp(-0.5*(AfMax-theParam.GetA1())*(AfMax-theParam.GetA1())/
00372                                (theParam.GetSigma1()*theParam.GetSigma1()));
00373 
00374     G4double P2 = std::exp(-0.5*(AfMax-theParam.GetA2())*(AfMax-theParam.GetA2())/
00375                            (theParam.GetSigma2()*theParam.GetSigma2()));
00376 
00377     Pas = P1+P2;
00378   }
00379 
00380   G4double Ps;
00381   if (theParam.GetW() < 0.001) Ps = 0.0;
00382   else {
00383     Ps = theParam.GetW()*std::exp(-0.5*(AfMax-theParam.GetAs())*(AfMax-theParam.GetAs())/
00384                                   (theParam.GetSigmaS()*theParam.GetSigmaS()));
00385   }
00386   G4double Psy = Ps/(Pas+Ps);
00387 
00388   // Fission fractions Xsy and Xas formed in symmetric and asymmetric modes
00389   G4double PPas = theParam.GetSigma1() + 2.0 * theParam.GetSigma2();
00390   G4double PPsy = theParam.GetW() * theParam.GetSigmaS();
00391   G4double Xas = PPas / (PPas+PPsy);
00392   G4double Xsy = PPsy / (PPas+PPsy);
00393 
00394   // Average kinetic energy for symmetric and asymmetric components
00395   G4double Eaverage = 0.1071*MeV*(Z*Z)/G4Pow::GetInstance()->Z13(A) + 22.2*MeV;
00396 
00397 
00398   // Compute maximal average kinetic energy of fragments and Energy Dispersion (sqrt)
00399   G4double TaverageAfMax;
00400   G4double ESigma;
00401   // Select randomly fission mode (symmetric or asymmetric)
00402   if (G4UniformRand() > Psy) { // Asymmetric Mode
00403     G4double A11 = theParam.GetA1()-0.7979*theParam.GetSigma1();
00404     G4double A12 = theParam.GetA1()+0.7979*theParam.GetSigma1();
00405     G4double A21 = theParam.GetA2()-0.7979*theParam.GetSigma2();
00406     G4double A22 = theParam.GetA2()+0.7979*theParam.GetSigma2();
00407     // scale factor
00408     G4double ScaleFactor = 0.5*theParam.GetSigma1()*(AsymmetricRatio(A,A11)+AsymmetricRatio(A,A12))+
00409       theParam.GetSigma2()*(AsymmetricRatio(A,A21)+AsymmetricRatio(A,A22));
00410     // Compute average kinetic energy for fragment with AfMax
00411     TaverageAfMax = (Eaverage + 12.5 * Xsy) * (PPas/ScaleFactor) * AsymmetricRatio(A,AfMax);
00412     ESigma = 10.0*MeV; // MeV
00413 
00414   } else { // Symmetric Mode
00415     G4double As0 = theParam.GetAs() + 0.7979*theParam.GetSigmaS();
00416     // scale factor
00417     G4double ScaleFactor = theParam.GetW()*theParam.GetSigmaS()*SymmetricRatio(A,As0);
00418     // Compute average kinetic energy for fragment with AfMax
00419     TaverageAfMax = (Eaverage - 12.5*MeV*Xas) * (PPsy/ScaleFactor) * SymmetricRatio(A,AfMax);
00420     ESigma = 8.0*MeV;
00421   }
00422 
00423 
00424   // Select randomly, in accordance with Gaussian distribution, fragment kinetic energy
00425   G4double KineticEnergy;
00426   G4int i = 0;
00427   do {
00428     KineticEnergy = G4RandGauss::shoot(TaverageAfMax,ESigma);
00429     if (i++ > 100) return Eaverage;
00430   } while (KineticEnergy < Eaverage-3.72*ESigma || 
00431            KineticEnergy > Eaverage+3.72*ESigma ||
00432            KineticEnergy > Tmax);
00433   
00434   return KineticEnergy;
00435 }
00436 
00437 G4double G4CompetitiveFission::AsymmetricRatio(G4int A, G4double A11)
00438 {
00439   const G4double B1 = 23.5;
00440   const G4double A00 = 134.0;
00441   return Ratio(G4double(A),A11,B1,A00);
00442 }
00443 
00444 G4double G4CompetitiveFission::SymmetricRatio(G4int A, G4double A11)
00445 {
00446   const G4double B1 = 5.32;
00447   const G4double A00 = A/2.0;
00448   return Ratio(G4double(A),A11,B1,A00);
00449 }
00450 
00451 G4double G4CompetitiveFission::Ratio(G4double A, G4double A11,
00452                                      G4double B1, G4double A00) 
00453 {
00454   if (A == 0.0) {
00455     throw G4HadronicException(__FILE__, __LINE__, 
00456                               "G4CompetitiveFission::Ratio: A == 0!");
00457   }
00458   if (A11 >= A/2.0 && A11 <= (A00+10.0)) {
00459     return 1.0-B1*((A11-A00)/A)*((A11-A00)/A);
00460   } else {
00461     return 1.0-B1*(10.0/A)*(10.0/A)-2.0*(10.0/A)*B1*((A11-A00-10.0)/A);
00462   }
00463 }
00464 
00465 G4ThreeVector G4CompetitiveFission::IsotropicVector(const G4double Magnitude)
00466   // Samples a isotropic random vectorwith a magnitud given by Magnitude.
00467   // By default Magnitude = 1.0
00468 {
00469   G4double CosTheta = 1.0 - 2.0*G4UniformRand();
00470   G4double SinTheta = std::sqrt(1.0 - CosTheta*CosTheta);
00471   G4double Phi = twopi*G4UniformRand();
00472   G4ThreeVector Vector(Magnitude*std::cos(Phi)*SinTheta,
00473                        Magnitude*std::sin(Phi)*SinTheta,
00474                        Magnitude*CosTheta);
00475   return Vector;
00476 }
00477 
00478 #ifdef debug
00479 void G4CompetitiveFission::CheckConservation(const G4Fragment & theInitialState,
00480                                              G4FragmentVector * Result) const
00481 {
00482     G4double ProductsEnergy =0;
00483     G4ThreeVector ProductsMomentum;
00484     G4int ProductsA = 0;
00485     G4int ProductsZ = 0;
00486     G4FragmentVector::iterator h;
00487     for (h = Result->begin(); h != Result->end(); h++) {
00488         G4LorentzVector tmp = (*h)->GetMomentum();
00489         ProductsEnergy += tmp.e();
00490         ProductsMomentum += tmp.vect();
00491         ProductsA += (*h)->GetA_asInt();
00492         ProductsZ += (*h)->GetZ_asInt();
00493     }
00494 
00495     if (ProductsA != theInitialState.GetA_asInt()) {
00496         G4cout << "!!!!!!!!!! Baryonic Number Conservation Violation !!!!!!!!!!" << G4endl;
00497         G4cout << "G4CompetitiveFission.cc: Barionic Number Conservation test for fission fragments" 
00498                << G4endl; 
00499         G4cout << "Initial A = " << theInitialState.GetA_asInt() 
00500                << "   Fragments A = " << ProductsA << "   Diference --> " 
00501                << theInitialState.GetA_asInt() - ProductsA << G4endl;
00502     }
00503     if (ProductsZ != theInitialState.GetZ_asInt()) {
00504         G4cout << "!!!!!!!!!! Charge Conservation Violation !!!!!!!!!!" << G4endl;
00505         G4cout << "G4CompetitiveFission.cc: Charge Conservation test for fission fragments" 
00506                << G4endl; 
00507         G4cout << "Initial Z = " << theInitialState.GetZ_asInt() 
00508                << "   Fragments Z = " << ProductsZ << "   Diference --> " 
00509                << theInitialState.GetZ() - ProductsZ << G4endl;
00510     }
00511     if (std::fabs(ProductsEnergy-theInitialState.GetMomentum().e()) > 1.0*keV) {
00512         G4cout << "!!!!!!!!!! Energy Conservation Violation !!!!!!!!!!" << G4endl;
00513         G4cout << "G4CompetitiveFission.cc: Energy Conservation test for fission fragments" 
00514                << G4endl; 
00515         G4cout << "Initial E = " << theInitialState.GetMomentum().e()/MeV << " MeV"
00516                << "   Fragments E = " << ProductsEnergy/MeV  << " MeV   Diference --> " 
00517                << (theInitialState.GetMomentum().e() - ProductsEnergy)/MeV << " MeV" << G4endl;
00518     } 
00519     if (std::fabs(ProductsMomentum.x()-theInitialState.GetMomentum().x()) > 1.0*keV || 
00520         std::fabs(ProductsMomentum.y()-theInitialState.GetMomentum().y()) > 1.0*keV ||
00521         std::fabs(ProductsMomentum.z()-theInitialState.GetMomentum().z()) > 1.0*keV) {
00522         G4cout << "!!!!!!!!!! Momentum Conservation Violation !!!!!!!!!!" << G4endl;
00523         G4cout << "G4CompetitiveFission.cc: Momentum Conservation test for fission fragments" 
00524                << G4endl; 
00525         G4cout << "Initial P = " << theInitialState.GetMomentum().vect() << " MeV"
00526                << "   Fragments P = " << ProductsMomentum  << " MeV   Diference --> " 
00527                << theInitialState.GetMomentum().vect() - ProductsMomentum << " MeV" << G4endl;
00528     }
00529     return;
00530 }
00531 #endif
00532 
00533 
00534 
00535 

Generated on Mon May 27 17:47:56 2013 for Geant4 by  doxygen 1.4.7