Changeset 1055 for trunk/source/processes/hadronic/models/de_excitation/fission/src/G4CompetitiveFission.cc

Timestamp:

May 28, 2009, 4:26:57 PM (15 years ago)

Author:

garnier

Message:

maj sur la beta de geant 4.9.3

File:

• trunk/source/processes/hadronic/models/de_excitation/fission/src/G4CompetitiveFission.cc (modified) (8 diffs)

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

@@ -25 +25 @@
 //
 //
-// $Id: G4CompetitiveFission.cc,v 1.9 2008/11/20 13:46:27 dennis Exp $
-// GEANT4 tag $Name: geant4-09-02 $
+// $Id: G4CompetitiveFission.cc,v 1.11 2009/03/13 18:57:17 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
 //
 // Hadronic Process: Nuclear De-excitations
 // by V. Lara (Oct 1998)
+//
+// J. M. Quesada (March 2009). Bugs fixed:
+//          - Full relativistic calculation (Lorentz boosts)
+//          - Fission pairing energy is included in fragment excitation energies
+// Now Energy and momentum are conserved in fission 
 
 #include "G4CompetitiveFission.hh"
@@ -127 +132 @@
     }
 
+
     // Atomic Mass of Nucleus (in MeV)
-    G4double M = G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(Z,A)/MeV;
+    G4double M = G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(Z,A);
+
     // Nucleus Momentum
     G4LorentzVector theNucleusMomentum = theNucleus.GetMomentum();
@@ -147 +154 @@
     G4double FragmentsExcitationEnergy = 0.0;
     G4double FragmentsKineticEnergy = 0.0;
+
+    //JMQ 04/03/09 It will be used latter to fix the bug in energy conservation
+    G4double FissionPairingEnergy=G4PairingCorrection::GetInstance()->GetFissionPairingCorrection(A,Z);
 
     G4int Trials = 0;
@@ -162 +172 @@
         if (A2 < 1 || Z2 < 0) 
             throw G4HadronicException(__FILE__, __LINE__, "G4CompetitiveFission::BreakUp: Can't define second fragment! ");
-        M2 = G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(Z2,A2)/MeV;
+
+        M2 = G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(Z2,A2);
 
         // Check that fragment masses are less or equal than total energy
-        //  if (M1 + M2 > theNucleusMomentum.mag()/MeV)
-        if (M1 + M2 > theNucleusMomentum.e()/MeV)
+        if (M1 + M2 > theNucleusMomentum.e())
             throw G4HadronicException(__FILE__, __LINE__, "G4CompetitiveFission::BreakUp: Fragments Mass > Total Energy");
 
         // Maximal Kinetic Energy (available energy for fragments)
-        //  G4double Tmax = theNucleusMomentum.mag()/MeV - M1 - M2;
         G4double Tmax = M + U - M1 - M2;
 
@@ -180 +189 @@
     
         // Excitation Energy
-        FragmentsExcitationEnergy = Tmax - FragmentsKineticEnergy;
+        //      FragmentsExcitationEnergy = Tmax - FragmentsKineticEnergy;
+        // JMQ 04/03/09 BUG FIXED: in order to fulfill energy conservation the
+        // fragments carry the fission pairing energy in form of 
+        //excitation energy
+
+        FragmentsExcitationEnergy = Tmax - FragmentsKineticEnergy+FissionPairingEnergy;
+
+    } while (FragmentsExcitationEnergy < 0.0 && Trials++ < 100);
     
-    } while (FragmentsExcitationEnergy < 0.0 && Trials++ < 100);
-  
   
 
@@ -190 +204 @@
   
 
-  // while (FragmentsExcitationEnergy < 0 && Trials < 100);
-  
-  // Fragment 1
+    // while (FragmentsExcitationEnergy < 0 && Trials < 100);
+  
+    // Fragment 1
     G4double U1 = FragmentsExcitationEnergy * (static_cast<G4double>(A1)/static_cast<G4double>(A));
     // Fragment 2
@@ -198 +212 @@
 
 
-    G4double Pmax = std::sqrt( 2 * ( ( (M1+U1)*(M2+U2) ) /
-                                ( (M1+U1)+(M2+U2) ) ) * FragmentsKineticEnergy);
-
-    G4ParticleMomentum momentum1 = IsotropicVector( Pmax );
-    G4ParticleMomentum momentum2( -momentum1 );
+    //JMQ  04/03/09 Full relativistic calculation is performed
+    //
+    G4double Fragment1KineticEnergy=(FragmentsKineticEnergy*(FragmentsKineticEnergy+2*(M2+U2)))/(2*(M1+U1+M2+U2+FragmentsKineticEnergy));
+    G4ParticleMomentum Momentum1(IsotropicVector(std::sqrt(Fragment1KineticEnergy*(Fragment1KineticEnergy+2*(M1+U1)))));
+    G4ParticleMomentum Momentum2(-Momentum1);
+    G4LorentzVector FourMomentum1(Momentum1,std::sqrt(Momentum1.mag2()+(M1+U1)*(M1+U1)));
+    G4LorentzVector FourMomentum2(Momentum2,std::sqrt(Momentum2.mag2()+(M2+U2)*(M2+U2)));
+
+    //JMQ 04/03/09 now we do Lorentz boosts (instead of Galileo boosts)
+    FourMomentum1.boost(theNucleusMomentum.boostVector());
+    FourMomentum2.boost(theNucleusMomentum.boostVector());
+    
+    //////////JMQ 04/03: Old version calculation is commented
+    // There was vioation of energy momentum conservation
+
+    //    G4double Pmax = std::sqrt( 2 * ( ( (M1+U1)*(M2+U2) ) /
+    //                          ( (M1+U1)+(M2+U2) ) ) * FragmentsKineticEnergy);
+
+    //G4ParticleMomentum momentum1 = IsotropicVector( Pmax );
+    //  G4ParticleMomentum momentum2( -momentum1 );
 
     // Perform a Galileo boost for fragments
-    momentum1 += (theNucleusMomentum.boostVector() * (M1+U1));
-    momentum2 += (theNucleusMomentum.boostVector() * (M2+U2));
+    //    momentum1 += (theNucleusMomentum.boostVector() * (M1+U1));
+    //    momentum2 += (theNucleusMomentum.boostVector() * (M2+U2));
 
 
     // Create 4-momentum for first fragment
     // Warning!! Energy conservation is broken
-    G4LorentzVector FourMomentum1( momentum1 , std::sqrt(momentum1.mag2() + (M1+U1)*(M1+U1)));
+    //JMQ 04/03/09 ...NOT ANY MORE!! BUGS FIXED: Energy and momentum are NOW conserved 
+    //    G4LorentzVector FourMomentum1( momentum1 , std::sqrt(momentum1.mag2() + (M1+U1)*(M1+U1)));
 
     // Create 4-momentum for second fragment
     // Warning!! Energy conservation is broken
-    G4LorentzVector FourMomentum2( momentum2 , std::sqrt(momentum2.mag2() + (M2+U2)*(M2+U2)));
+    //JMQ 04/03/09 ...NOT ANY MORE!! BUGS FIXED: Energy and momentum are NOW conserved
+    //    G4LorentzVector FourMomentum2( momentum2 , std::sqrt(momentum2.mag2() + (M2+U2)*(M2+U2)));
+
+    //////////
 
     // Create Fragments
@@ -227 +260 @@
     Fragment2->SetCreatorModel(G4String("G4CompetitiveFission"));
 #endif
-  // Create Fragment Vector
+    // Create Fragment Vector
     G4FragmentVector * theResult = new G4FragmentVector;