Changeset 1340 for trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/src

Timestamp:

Nov 5, 2010, 3:45:55 PM (14 years ago)

Author:

garnier

Message:

update ti head

Location:

trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/src

Files:

• G4GNASHTransitions.cc (modified) (3 diffs)
• G4HETCAlpha.cc (modified) (2 diffs)
• G4HETCChargedFragment.cc (modified) (1 diff)
• G4HETCDeuteron.cc (modified) (2 diffs)
• G4HETCEmissionFactory.cc (modified) (2 diffs)
• G4HETCFragment.cc (modified) (6 diffs)
• G4HETCHe3.cc (modified) (2 diffs)
• G4HETCNeutron.cc (modified) (2 diffs)
• G4HETCProton.cc (modified) (1 diff)
• G4HETCTriton.cc (modified) (2 diffs)
• G4LowEIonFragmentation.cc (modified) (6 diffs)
• G4PreCompoundAlpha.cc (modified) (12 diffs)
• G4PreCompoundDeuteron.cc (modified) (12 diffs)
• G4PreCompoundEmission.cc (modified) (14 diffs)
• G4PreCompoundEmissionFactory.cc (modified) (2 diffs)
• G4PreCompoundFragment.cc (modified) (2 diffs)
• G4PreCompoundFragmentVector.cc (modified) (1 diff)
• G4PreCompoundHe3.cc (modified) (11 diffs)
• G4PreCompoundIon.cc (modified) (3 diffs)
• G4PreCompoundModel.cc (modified) (7 diffs)
• G4PreCompoundNeutron.cc (modified) (11 diffs)
• G4PreCompoundNucleon.cc (modified) (2 diffs)
• G4PreCompoundParameters.cc (modified) (1 diff)
• G4PreCompoundProton.cc (modified) (14 diffs)
• G4PreCompoundTransitions.cc (modified) (6 diffs)
• G4PreCompoundTriton.cc (modified) (10 diffs)
• G4VPreCompoundEmissionFactory.cc (modified) (1 diff)
• G4VPreCompoundFragment.cc (modified) (5 diffs)

trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/src/G4GNASHTransitions.cc

@@ -24 +24 @@
 // ********************************************************************
 //
+// $Id: G4GNASHTransitions.cc,v 1.6 2010/08/20 07:42:19 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03-ref-09 $
+//
+// 20.08.2010 V.Ivanchenko move constructor and destructor to the source 
+
 #include "G4GNASHTransitions.hh"
 #include "G4PreCompoundParameters.hh"
 #include "G4HadronicException.hh"
 
-const G4GNASHTransitions & G4GNASHTransitions::
-operator=(const G4GNASHTransitions & )
-{
-  throw G4HadronicException(__FILE__, __LINE__, "G4GNASHTransitions::operator= meant to not be accesable");
-  return *this;
-}
+G4GNASHTransitions::G4GNASHTransitions()
+{}
 
-G4bool G4GNASHTransitions::operator==(const G4GNASHTransitions & ) const
-{
-  return false;
-}
-
-G4bool G4GNASHTransitions::operator!=(const G4GNASHTransitions & ) const
-{
-  return true;
-}
-
+G4GNASHTransitions::~G4GNASHTransitions()
+{}
 
 G4double G4GNASHTransitions::
@@ -72 +65 @@
   Probability *= theMatrixElement;
 
-
   return Probability;
 }
 
-G4Fragment G4GNASHTransitions::
-PerformTransition(const G4Fragment & aFragment)
+void G4GNASHTransitions::PerformTransition(G4Fragment & result)
 {
-  G4Fragment result(aFragment);
   result.SetNumberOfParticles(result.GetNumberOfParticles()+1);
   result.SetNumberOfHoles(result.GetNumberOfHoles()+1);
@@ -91 +81 @@
       result.SetNumberOfCharged(result.GetNumberOfParticles());
     }
-
-  return result;
 }

trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/src/G4HETCAlpha.cc

@@ -24 +24 @@
 // ********************************************************************
 //
+// $Id: G4HETCAlpha.cc,v 1.4 2010/08/28 15:16:55 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03-ref-09 $
+//
+// by V. Lara
+//
+// Modified:
+// 23.08.2010 V.Ivanchenko general cleanup, move constructor and destructor 
+//            the source, use G4Pow
+
 #include "G4HETCAlpha.hh"
+#include "G4Alpha.hh"
+
+G4HETCAlpha::G4HETCAlpha() 
+  : G4HETCChargedFragment(G4Alpha::Alpha(), &theAlphaCoulombBarrier)
+{}
+
+G4HETCAlpha::~G4HETCAlpha() 
+{}
+
+G4double G4HETCAlpha::GetAlpha()
+{
+  G4double C = 0.0;
+  G4int aZ = GetZ() + GetRestZ();
+  if (aZ <= 30) 
+    {
+      C = 0.10;
+    } 
+  else if (aZ <= 50) 
+    {
+      C = 0.1 + -((aZ-50.)/20.)*0.02;
+    }
+  else if (aZ < 70) 
+    {
+      C = 0.08 + -((aZ-70.)/20.)*0.02;
+    }
+  else 
+    {
+      C = 0.06;
+    }
+  return 1.0+C;
+}
+  
+G4double G4HETCAlpha::GetBeta()
+{
+  return -GetCoulombBarrier();
+}
+  
+G4double G4HETCAlpha::GetSpinFactor()
+{
+  return 1.0;
+}
 
 G4double G4HETCAlpha::K(const G4Fragment & aFragment)
 {
-  if (GetStage() != 1) return 1.0;
-  // Number of protons in projectile
-  G4double Pa = static_cast<G4int>(aFragment.GetParticleDefinition()->GetPDGCharge());
-  // Number of neutrons in projectile 
-  G4double Na = aFragment.GetParticleDefinition()->GetBaryonNumber();
-  G4double TargetA = aFragment.GetA() - Na;
-  G4double TargetZ = aFragment.GetZ() - Pa;
-  Na -= Pa;
-  G4double r = TargetZ/TargetA;
+  // Number of protons in emitted fragment
+  G4int Pa = GetZ();
+  // Number of neutrons in emitted fragment 
+  G4int Na = GetA() - Pa;
 
-  
-  G4double P = aFragment.GetNumberOfParticles();
-  G4double H = aFragment.GetNumberOfHoles();
+  G4int TargetZ = GetRestZ();
+  G4int TargetA = GetRestA();
+  G4double r = G4double(TargetZ)/G4double(TargetA);
+
+  G4int P = aFragment.GetNumberOfParticles();
+  G4int H = aFragment.GetNumberOfHoles();
 
   G4double result = 0.0;
@@ -52 +100 @@
          Pa*(Pa-1.0)*Na*(Na-1.0));
 
-      result /= 6.0*std::pow((TargetZ/TargetA)*((TargetA-TargetZ)/TargetA),2.0);
+      result /= 6.0*r*r*(1. - r) *(1. - r);
     }
-
   return std::max(0.0,result);
-
 }

trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/src/G4HETCChargedFragment.cc

@@ -24 +24 @@
 // ********************************************************************
 //
+// $Id: G4HETCChargedFragment.cc,v 1.3 2010/08/28 15:16:55 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03-ref-09 $
+//
 // by V. Lara
+//
+// Modified:
+// 23.08.2010 V.Ivanchenko general cleanup, move constructor and destructor 
+//            the source, use G4Pow
+//
 
 #include "G4HETCChargedFragment.hh"
+#include "G4VCoulombBarrier.hh"
 #include "G4PreCompoundParameters.hh"
 
+G4HETCChargedFragment::G4HETCChargedFragment(
+  const G4ParticleDefinition* pd, G4VCoulombBarrier * aCoulombBarrier)
+  : G4HETCFragment(pd, aCoulombBarrier)
+{}
+
+G4HETCChargedFragment::~G4HETCChargedFragment()
+{}
 
 G4double G4HETCChargedFragment::
 GetKineticEnergy(const G4Fragment & aFragment)
 {
-  // Number of protons in projectile
-  G4double Pa = aFragment.GetParticleDefinition()->GetPDGCharge();
-  // Number of neutrons in projectile 
-  G4double Na = aFragment.GetParticleDefinition()->GetBaryonNumber();
-  Na -= Pa;
-  
-  G4double Pb = aFragment.GetNumberOfParticles();
-  G4double H = aFragment.GetNumberOfHoles();
+  G4int Pb = aFragment.GetNumberOfParticles();
+  G4int H = aFragment.GetNumberOfHoles();
 
-  G4double Ab = std::max(0.0,(Pb*Pb+H*H+Pb-3*H)/4.0);
+  G4double g0 = (6.0/pi2)*aFragment.GetA_asInt()*theParameters->GetLevelDensity();
+
+  G4double Ab = std::max(0.0,G4double(Pb*Pb+H*H+Pb-3*H)/(4.0*g0));
   G4double Emax = GetMaximalKineticEnergy() - Ab;
 
-  G4double x = BetaRand(static_cast<G4int>(Pb+H),2);
+  G4double x = BetaRand(Pb + H, 2);
   
   return Emax - (Emax-GetCoulombBarrier())*x;

trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/src/G4HETCDeuteron.cc

@@ -24 +24 @@
 // ********************************************************************
 //
+// $Id: G4HETCDeuteron.cc,v 1.3 2010/08/28 15:16:55 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03-ref-09 $
+//
+// by V. Lara
+//
+// Modified:
+// 23.08.2010 V.Ivanchenko general cleanup, move constructor and destructor 
+//            the source, use G4Pow
+
 #include "G4HETCDeuteron.hh"
+#include "G4Deuteron.hh"
+
+G4HETCDeuteron::G4HETCDeuteron() 
+  : G4HETCChargedFragment(G4Deuteron::Deuteron(), &theDeuteronCoulombBarrier)
+{}
+
+G4HETCDeuteron::~G4HETCDeuteron() 
+{}
+
+G4double G4HETCDeuteron::GetAlpha()
+{
+  G4double C = 0.0;
+  G4int aZ = GetZ() + GetRestZ();
+  if (aZ >= 70) 
+    {
+      C = 0.10;
+    } 
+  else 
+    {
+      C = ((((0.15417e-06*aZ) - 0.29875e-04)*aZ + 0.21071e-02)*aZ - 0.66612e-01)*aZ + 0.98375; 
+    }
+  return 1.0 + C/2.0;
+}
+  
+G4double G4HETCDeuteron::GetBeta()
+{
+  return -GetCoulombBarrier();
+}
+
+G4double G4HETCDeuteron::GetSpinFactor()
+{
+  // 2s+1
+  return 3.0;
+}
 
 G4double G4HETCDeuteron::K(const G4Fragment & aFragment)
 {
-  if (GetStage() != 1) return 1.0;
-  // Number of protons in projectile
-  G4double Pa = static_cast<G4int>(aFragment.GetParticleDefinition()->GetPDGCharge());
-  // Number of neutrons in projectile 
-  G4double Na = aFragment.GetParticleDefinition()->GetBaryonNumber();
-  G4double TargetA = aFragment.GetA() - Na;
-  G4double TargetZ = aFragment.GetZ() - Pa;
-  Na -= Pa;
-  G4double r = TargetZ/TargetA;
+  // Number of protons in emitted fragment
+  G4int Pa = GetZ();
+  // Number of neutrons in emitted fragment 
+  G4int Na = GetA() - Pa;
 
+  G4int TargetZ = GetRestZ();
+  G4int TargetA = GetRestA();
+  G4double r = G4double(TargetZ)/G4double(TargetA);
   
-  G4double P = aFragment.GetNumberOfParticles();
-  G4double H = aFragment.GetNumberOfHoles();
+  G4int P = aFragment.GetNumberOfParticles();
+  G4int H = aFragment.GetNumberOfHoles();
 
   G4double result = 0.0;
@@ -47 +88 @@
       result = 2.0* (H*(H-1.0)*r*(r-1.0)+H*(Na*r+Pa*(1.0-r)) + Pa*Na)/(P*(P-1.0));
 
-      result /= (TargetZ/TargetA)*((TargetA-TargetZ)/TargetA);
+      result /= r*(1.0 - r);
     }
-
   return std::max(0.0,result);
-
 }

trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/src/G4HETCEmissionFactory.cc

@@ -24 +24 @@
 // ********************************************************************
 //
+// $Id: G4HETCEmissionFactory.cc,v 1.5 2010/08/28 15:16:55 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03-ref-09 $
+//
+// by V. Lara
+//
+// Modified:
+// 23.08.2010 V.Ivanchenko general cleanup, move constructor and destructor 
+//            the source
+//
+
 #include "G4HETCEmissionFactory.hh"
 
@@ -33 +43 @@
 #include "G4HETCAlpha.hh"
 
+G4HETCEmissionFactory::G4HETCEmissionFactory()
+{}
 
-const G4HETCEmissionFactory & G4HETCEmissionFactory::
-operator=(const G4HETCEmissionFactory & )
-{
-  throw G4HadronicException(__FILE__, __LINE__, "G4HETCEmissionFactory::operator= meant to not be accessable.");
-  return *this;
-}
-
-G4bool G4HETCEmissionFactory::
-operator==(const G4HETCEmissionFactory & ) const
-{
-  throw G4HadronicException(__FILE__, __LINE__, "G4HETCEmissionFactory::operator== meant to not be accessable.");
-  return false;
-}
-
-G4bool G4HETCEmissionFactory::
-operator!=(const G4HETCEmissionFactory & ) const
-{
-  throw G4HadronicException(__FILE__, __LINE__, "G4HETCEmissionFactory::operator!= meant to not be accessable.");
-  return true;
-}
-
+G4HETCEmissionFactory::~G4HETCEmissionFactory()
+{}
 
 std::vector<G4VPreCompoundFragment*> *  G4HETCEmissionFactory::

trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/src/G4HETCFragment.cc

@@ -24 +24 @@
 // ********************************************************************
 //
+// $Id: G4HETCFragment.cc,v 1.4 2010/08/28 15:16:55 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03-ref-09 $
+//
 // by V. Lara
+//
+// Modified:
+// 23.08.2010 V.Ivanchenko general cleanup, move constructor and destructor 
+//            the source, use G4Pow
  
 #include "G4HETCFragment.hh"
@@ -30 +37 @@
 
 G4HETCFragment::
-G4HETCFragment(const G4HETCFragment & right) :
-  G4VPreCompoundFragment(right)
+G4HETCFragment(const G4ParticleDefinition* part,
+               G4VCoulombBarrier* aCoulombBarrier)
+  : G4VPreCompoundFragment(part, aCoulombBarrier)
 {
+  G4double r0 = theParameters->Getr0();
+  r2norm = r0*r0/(CLHEP::pi*CLHEP::hbarc*CLHEP::hbarc*CLHEP::hbarc);
 }
 
-G4HETCFragment::
-G4HETCFragment(const G4double anA,
-               const G4double aZ, 
-               G4VCoulombBarrier* aCoulombBarrier,
-               const G4String & aName) :
-  G4VPreCompoundFragment(anA,aZ,aCoulombBarrier,aName)
+G4HETCFragment::~G4HETCFragment()
 {}
-
-G4HETCFragment::~G4HETCFragment()
-{
-}
-
-const G4HETCFragment & G4HETCFragment::
-operator= (const G4HETCFragment & right)
-{
-  if (&right != this) this->G4VPreCompoundFragment::operator=(right);
-  return *this;
-}
-
-G4int G4HETCFragment::operator==(const G4HETCFragment & right) const
-{
-  return G4VPreCompoundFragment::operator==(right);
-}
-
-G4int G4HETCFragment::operator!=(const G4HETCFragment & right) const
-{
-  return G4VPreCompoundFragment::operator!=(right);
-}
-
-
 
 G4double G4HETCFragment::
@@ -70 +52 @@
 {
   if (GetEnergyThreshold() <= 0.0) 
-  {
+    {
       theEmissionProbability = 0.0;
       return 0.0;
-  }    
+    }    
   // Coulomb barrier is the lower limit 
   // of integration over kinetic energy
@@ -80 +62 @@
   // Excitation energy of nucleus after fragment emission is the upper limit
   // of integration over kinetic energy
-  G4double UpperLimit = this->GetMaximalKineticEnergy();
+  G4double UpperLimit = GetMaximalKineticEnergy();
   
-  theEmissionProbability = IntegrateEmissionProbability(LowerLimit,UpperLimit,aFragment);
+  theEmissionProbability = 
+    IntegrateEmissionProbability(LowerLimit,UpperLimit,aFragment);
     
   return theEmissionProbability;
@@ -92 +75 @@
 {
   
-  if ( !IsItPossible(aFragment) ) return 0.0;
-
-  const G4double r0 = G4PreCompoundParameters::GetAddress()->Getr0();
+  if ( !IsItPossible(aFragment) ) { return 0.0; }
     
   G4double U = aFragment.GetExcitationEnergy();
-  G4double P = aFragment.GetNumberOfParticles();
-  G4double H = aFragment.GetNumberOfHoles();
-  G4double N = P + H;
-  G4double Pb = P - GetA();
-  G4double Nb = Pb + H;
-  if (Nb <= 0.0) return 0.0;
-  
-  G4double A = (P*P+H*H+P-3*H)/4.0;
-  G4double Ab = (Pb*Pb+H*H+Pb-3*H)/4.0;
+
+  G4int P  = aFragment.GetNumberOfParticles();
+  G4int H  = aFragment.GetNumberOfHoles();
+  G4int N  = P + H;
+  G4int Pb = P - GetA();
+  G4int Nb = Pb + H;
+  if (Nb <= 0.0) { return 0.0; }
+  G4double g = (6.0/pi2)*aFragment.GetA()*theParameters->GetLevelDensity();
+  G4double gb = (6.0/pi2)*GetRestA()*theParameters->GetLevelDensity();
+
+  G4double A  = G4double(P*P+H*H+P-3*H)/(4.0*g);
+  G4double Ab = G4double(Pb*Pb+H*H+Pb-3*H)/(4.0*gb);
   U = std::max(U-A,0.0);
-  if (U <= 0.0) return 0.0;
+  if (U <= 0.0) { return 0.0; }
 
-  G4double g = (6.0/pi2)*aFragment.GetA()*G4PreCompoundParameters::GetAddress()->GetLevelDensity();
-  G4double gb = (6.0/pi2)*GetRestA()*G4PreCompoundParameters::GetAddress()->GetLevelDensity();
-
-  G4double Pf = P;
-  G4double Hf = H;
-  G4double Nf = N-1.0;
-  for (G4int i = 1; i < GetA(); i++)
+  G4int Pf = P;
+  G4int Hf = H;
+  G4int Nf = N-1;
+  for (G4int i = 1; i < GetA(); ++i)
     {
       Pf *= (P-i);
@@ -125 +106 @@
   G4double Y = std::max(Up - Ab - Low, 0.0);
 
-  G4double Probability = GetSpinFactor()/(pi*hbarc*hbarc*hbarc) * GetReducedMass() * GetAlpha() *
-    r0 * r0 * std::pow(GetRestA(),2.0/3.0)/std::pow(U,N-1) * (std::pow(gb,Nb)/std::pow(g,N)) * Pf * Hf * Nf * K(aFragment) *
-    std::pow(Y,Nb) * (X/Nb - Y/(Nb+1));
+  G4double Probability = r2norm*GetSpinFactor()*GetReducedMass()*GetAlpha() 
+    *g4pow->Z23(GetRestA())*Pf*Hf*Nf*K(aFragment)*(X/Nb - Y/(Nb+1))
+    *U*g4pow->powN(g*gb,Nb)/g4pow->powN(g*U,N);
+
+  //  G4double Probability = GetSpinFactor()/(pi*hbarc*hbarc*hbarc) 
+  //  * GetReducedMass() * GetAlpha() *
+  //  r0 * r0 * std::pow->Z23(GetRestA())/std::pow->pow(U,G4double(N-1)) * 
+  //  (std::pow->(gb,Nb)/std::pow(g,N)) * Pf * Hf * Nf * K(aFragment) *
+  //  std::pow(Y,Nb) * (X/Nb - Y/(Nb+1));
 
   return Probability;

trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/src/G4HETCHe3.cc

@@ -24 +24 @@
 // ********************************************************************
 //
+// $Id: G4HETCHe3.cc,v 1.4 2010/08/28 15:16:55 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03-ref-09 $
+//
+// by V. Lara
+//
+// Modified:
+// 23.08.2010 V.Ivanchenko general cleanup, move constructor and destructor 
+//            the source, use G4Pow
+
 #include "G4HETCHe3.hh"
+#include "G4He3.hh"
+
+G4HETCHe3::G4HETCHe3() 
+  : G4HETCChargedFragment(G4He3::He3(), &theHe3CoulombBarrier)
+{}
+
+G4HETCHe3::~G4HETCHe3() 
+{}
+
+G4double G4HETCHe3::GetAlpha()
+{
+  G4double C = 0.0;
+  G4int aZ = GetZ() + GetRestZ();
+  if (aZ <= 30) 
+    {
+      C = 0.10;
+    } 
+  else if (aZ <= 50) 
+    {
+      C = 0.1 + -((aZ-50.)/20.)*0.02;
+    } 
+  else if (aZ < 70) 
+    {
+      C = 0.08 + -((aZ-70.)/20.)*0.02;
+    }
+  else 
+    {
+      C = 0.06;
+    }
+  return 1.0 + C*(4.0/3.0);
+}
+  
+G4double G4HETCHe3::GetBeta()
+{
+  return -GetCoulombBarrier();
+}
+
+G4double G4HETCHe3::GetSpinFactor()
+{
+  // 2s+1
+  return 2.0;
+}    
 
 G4double G4HETCHe3::K(const G4Fragment & aFragment)
 {
-  if (GetStage() != 1) return 1.0;
-  // Number of protons in projectile
-  G4double Pa = static_cast<G4int>(aFragment.GetParticleDefinition()->GetPDGCharge());
-  // Number of neutrons in projectile 
-  G4double Na = aFragment.GetParticleDefinition()->GetBaryonNumber();
-  G4double TargetA = aFragment.GetA() - Na;
-  G4double TargetZ = aFragment.GetZ() - Pa;
-  Na -= Pa;
-  G4double r = TargetZ/TargetA;
+  // Number of protons in emitted fragment
+  G4int Pa = GetZ();
+  // Number of neutrons in emitted fragment 
+  G4int Na = GetA() - Pa;
 
+  G4int TargetZ = GetRestZ();
+  G4int TargetA = GetRestA();
+  G4double r = G4double(TargetZ)/G4double(TargetA);
 
-  
-  G4double P = aFragment.GetNumberOfParticles();
-  G4double H = aFragment.GetNumberOfHoles();
+  G4int P = aFragment.GetNumberOfParticles();
+  G4int H = aFragment.GetNumberOfHoles();
 
   G4double result = 0.0;
@@ -52 +100 @@
          Pa*Na*(Pa-1.0));
 
-      result /= 3.0*std::pow(TargetZ/TargetA, 2.0) * ((TargetA-TargetZ)/TargetA);
+      result /= 3.0*r*r*(1.0 - r);
     }
-
   return std::max(0.0,result);
-
 }

trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/src/G4HETCNeutron.cc

@@ -24 +24 @@
 // ********************************************************************
 //
+// $Id: G4HETCNeutron.cc,v 1.3 2010/08/28 15:16:55 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03-ref-09 $
+//
+// by V. Lara
+//
+// Modified:
+// 23.08.2010 V.Ivanchenko general cleanup, move constructor and destructor 
+//            the source, use G4Pow
+
 #include "G4HETCNeutron.hh"
+#include "G4Neutron.hh"
 
+G4HETCNeutron::G4HETCNeutron() 
+  : G4HETCFragment(G4Neutron::Neutron(), &theNeutronCoulombBarrier)
+{}
+
+G4HETCNeutron::~G4HETCNeutron() 
+{}
+
+G4double G4HETCNeutron::GetAlpha()
+{
+  return 0.76+2.2/g4pow->Z13(GetRestA());
+}
+  
+G4double G4HETCNeutron::GetBeta() 
+{
+  return (2.12/g4pow->Z23(GetRestA())-0.05)*MeV/GetAlpha();
+}
+
+G4double G4HETCNeutron::GetSpinFactor()
+{
+  // (2s+1)
+  return 2.0;
+}
+  
 G4double G4HETCNeutron::K(const G4Fragment & aFragment)
 {
-  if (GetStage() != 1) return 1.0;
-  // Number of protons in projectile
-  G4double Pa = static_cast<G4int>(aFragment.GetParticleDefinition()->GetPDGCharge());
-  // Number of neutrons in projectile 
-  G4double Na = aFragment.GetParticleDefinition()->GetBaryonNumber();
-  G4double TargetA = aFragment.GetA() - Na;
-  G4double TargetZ = aFragment.GetZ() - Pa;
-  Na -= Pa;
-  G4double r = TargetZ/TargetA;
+  // Number of protons in emitted fragment
+  G4int Pa = GetZ();
+  // Number of neutrons in emitted fragment 
+  G4int Na = GetA() - Pa;
+
+  G4int TargetZ = GetRestZ();
+  G4int TargetA = GetRestA();
+  G4double r = G4double(TargetZ)/G4double(TargetA);
   
-  G4double P = aFragment.GetNumberOfParticles();
-  G4double H = aFragment.GetNumberOfHoles();
+  G4int P = aFragment.GetNumberOfParticles();
+  G4int H = aFragment.GetNumberOfHoles();
   
   G4double result = 0.0;
   if (P > 0)
     {
-      result = (H*(1.0-r)+Na)/P;
-
-      result /= (TargetA-TargetZ)/TargetA;
+      result = (H + Na/(1.0-r))/P;
     }
   
@@ -52 +82 @@
 }
 
-
 G4double G4HETCNeutron::GetKineticEnergy(const G4Fragment & aFragment)
 {
-  G4double H = aFragment.GetNumberOfHoles();
-  G4double Pb = aFragment.GetNumberOfParticles() - GetA();
-  G4double Nb = Pb + H;
+  G4int H = aFragment.GetNumberOfHoles();
+  G4int Pb = aFragment.GetNumberOfParticles();
+  G4int Nb = Pb + H;
+  G4double g0 = (6.0/pi2)*aFragment.GetA_asInt()*theParameters->GetLevelDensity();
   
-  G4double Ab = std::max(0.0,(Pb*Pb+H*H+Pb-3*H)/4.0);
+  G4double Ab = std::max(0.0,G4double(Pb*Pb+H*H+Pb-3*H)/(4.0*g0));
   G4double Emax = GetMaximalKineticEnergy() - Ab;
   
-  G4double cut = GetBeta() / (GetBeta()+Emax/(Nb+1));
+  G4double cut = GetBeta() / (GetBeta()+Emax/G4double(Nb+1));
   G4double x(0.0);
   if (G4UniformRand() <= cut)
     {
-      x = BetaRand(static_cast<G4int>(Nb),1);
+      x = BetaRand(Nb,1);
     }
   else 
     {
-      x = BetaRand(static_cast<G4int>(Nb),2);
+      x = BetaRand(Nb,2);
     }
 

trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/src/G4HETCProton.cc

@@ -24 +24 @@
 // ********************************************************************
 //
+// $Id: G4HETCProton.cc,v 1.3 2010/08/28 15:16:55 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03-ref-09 $
+//
+// by V. Lara
+//
+// Modified:
+// 23.08.2010 V.Ivanchenko general cleanup, move constructor and destructor 
+//            the source, use G4Pow
+
 #include "G4HETCProton.hh"
+#include "G4Proton.hh"
+
+G4HETCProton::G4HETCProton() 
+  : G4HETCChargedFragment(G4Proton::Proton(), &theProtonCoulombBarrier)
+{}
+
+G4HETCProton::~G4HETCProton() 
+{}
+
+G4double G4HETCProton::GetAlpha()
+{
+  G4int aZ = GetRestZ();
+  G4double C = 0.0;
+  if (aZ >= 70) 
+    {
+      C = 0.10;
+    } 
+  else 
+    {
+      C = ((((0.15417e-06*aZ) - 0.29875e-04)*aZ + 0.21071e-02)*aZ - 0.66612e-01)*aZ + 0.98375;
+    }
+  return 1.0 + C;
+}
+  
+G4double G4HETCProton::GetBeta()
+{
+  return -GetCoulombBarrier();
+}
+  
+G4double G4HETCProton::GetSpinFactor()
+{
+  // 2s+1
+  return 2.0;
+}
 
 G4double G4HETCProton::K(const G4Fragment & aFragment)
 {
-  if (GetStage() != 1) return 1.0;
-  // Number of protons in projectile
-  G4double Pa = static_cast<G4int>(aFragment.GetParticleDefinition()->GetPDGCharge());
-  // Number of neutrons in projectile 
-  G4double Na = aFragment.GetParticleDefinition()->GetBaryonNumber();
-  G4double TargetA = aFragment.GetA() - Na;
-  G4double TargetZ = aFragment.GetZ() - Pa;
-  Na -= Pa;
-  G4double r = TargetZ/TargetA;
+  // Number of protons in emitted fragment
+  G4int Pa = GetZ();
 
+  G4int TargetZ = GetRestZ();
+  G4int TargetA = GetRestA();
+  G4double r = G4double(TargetZ)/G4double(TargetA);
 
-  G4double P = aFragment.GetNumberOfParticles();
-  G4double H = aFragment.GetNumberOfHoles();
+  G4int P = aFragment.GetNumberOfParticles();
+  G4int H = aFragment.GetNumberOfHoles();
 
   G4double result = 0.0;
   if (P > 0)
     {
-      result = (H*r + Pa)/P;
-      
-      result /= TargetZ/TargetA;
+      result = (H*r + Pa)/P/r;
     }
 

trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/src/G4HETCTriton.cc

@@ -24 +24 @@
 // ********************************************************************
 //
+// $Id: G4HETCTriton.cc,v 1.4 2010/08/28 15:16:55 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03-ref-09 $
+//
+// by V. Lara
+//
+// Modified:
+// 23.08.2010 V.Ivanchenko general cleanup, move constructor and destructor 
+//            the source, use G4Pow
+
 #include "G4HETCTriton.hh"
+#include "G4Triton.hh"
+
+G4HETCTriton::G4HETCTriton() 
+  : G4HETCChargedFragment(G4Triton::Triton(), &theTritonCoulombBarrier)
+{}
+
+G4HETCTriton::~G4HETCTriton() 
+{}
+
+G4double G4HETCTriton::GetAlpha()
+{
+  G4double C = 0.0;
+  G4int aZ = GetZ() + GetRestZ();
+  if (aZ >= 70) 
+    {
+      C = 0.10;
+    } 
+  else 
+    {
+      C = ((((0.15417e-06*aZ) - 0.29875e-04)*aZ + 0.21071e-02)*aZ - 0.66612e-01)*aZ + 0.98375; 
+    }
+  
+  return 1.0 + C/3.0;
+}
+  
+G4double G4HETCTriton::GetBeta()
+{
+  return -GetCoulombBarrier();
+}
+
+G4double G4HETCTriton::GetSpinFactor()
+{
+  // 2s+1
+  return 2.0;
+}
 
 G4double G4HETCTriton::K(const G4Fragment & aFragment)
 {
-  if (GetStage() != 1) return 1.0;
-  // Number of protons in projectile
-  G4double Pa = static_cast<G4int>(aFragment.GetParticleDefinition()->GetPDGCharge());
-  // Number of neutrons in projectile 
-  G4double Na = aFragment.GetParticleDefinition()->GetBaryonNumber();
-  G4double TargetA = aFragment.GetA() - Na;
-  G4double TargetZ = aFragment.GetZ() - Pa;
-  Na -= Pa;
-  G4double r = TargetZ/TargetA;
+  // Number of protons in emitted fragment
+  G4int Pa = GetZ();
+  // Number of neutrons in emitted fragment 
+  G4int Na = GetA() - Pa;
 
-  
-  G4double P = aFragment.GetNumberOfParticles();
-  G4double H = aFragment.GetNumberOfHoles();
+  G4int TargetZ = GetRestZ();
+  G4int TargetA = GetRestA();
+  G4double r = G4double(TargetZ)/G4double(TargetA);
+
+  G4int P = aFragment.GetNumberOfParticles();
+  G4int H = aFragment.GetNumberOfHoles();
 
   G4double result = 0.0;
@@ -51 +93 @@
          Pa*Na*(Na-1.0));
 
-      result /= 3.0*(TargetZ/TargetA)*std::pow((TargetA-TargetZ)/TargetA,2.0);
+      result /= 3.0*r*(1.0 - r)*(1.0 - r);
     }
-
   return std::max(0.0,result);
-
 }

trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/src/G4LowEIonFragmentation.cc

@@ -24 +24 @@
 // ********************************************************************
 //
-//
-// $Id: G4LowEIonFragmentation.cc,v 1.5 2010/06/01 16:51:11 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-04-beta-01 $
+// $Id: G4LowEIonFragmentation.cc,v 1.7 2010/09/08 16:38:09 gunter Exp $
+// GEANT4 tag $Name: geant4-09-03-ref-09 $
 //
 //---------------------------------------------------------------------------
@@ -72 +71 @@
 
   // Get Target A, Z
-  G4double aTargetA = theNucleus.GetN();
-  G4double aTargetZ = theNucleus.GetZ();
+  G4int aTargetA = theNucleus.GetA_asInt();
+  G4int aTargetZ = theNucleus.GetZ_asInt();
 
   // Get Projectile A, Z
-  G4double aProjectileA = thePrimary.GetDefinition()->GetBaryonNumber();
-  G4double aProjectileZ = thePrimary.GetDefinition()->GetPDGCharge();
+  G4int aProjectileA = thePrimary.GetDefinition()->GetBaryonNumber();
+  G4int aProjectileZ = G4lrint(thePrimary.GetDefinition()->GetPDGCharge());
 
   // Get Maximum radius of both
@@ -131 +130 @@
     }
   }
-  hits ++;
+  hits++;
 
   // From target:
@@ -157 +156 @@
   G4double compoundEnergy = thePrimary.GetTotalEnergy()*particlesFromProjectile/aProjectileA;  
   G4double targetMass = G4ParticleTable::GetParticleTable()
-                        ->GetIonTable()->GetIonMass(static_cast<G4int>(aTargetZ) ,static_cast<G4int>(aTargetA));
+                        ->GetIonTable()->GetIonMass(aTargetZ, aTargetA);
   compoundEnergy += targetMass;
   G4LorentzVector fragment4Momentum(exciton3Momentum, compoundEnergy);
@@ -163 +162 @@
   // take the nucleons and fill the Fragments
   G4Fragment anInitialState;
-  anInitialState.SetA(aTargetA+particlesFromProjectile);
-  anInitialState.SetZ(aTargetZ+chargedFromProjectile);
+  anInitialState.SetZandA_asInt(aTargetZ+chargedFromProjectile,
+                                aTargetA+particlesFromProjectile);
   // M.A. Cortes fix
   //anInitialState.SetNumberOfParticles(particlesFromProjectile);
@@ -187 +186 @@
  
     G4Fragment initialState2;
-    initialState2.SetA(aProjectileA-particlesFromProjectile);
-    initialState2.SetZ(aProjectileZ-chargedFromProjectile);
+    initialState2.SetZandA_asInt(aProjectileZ-chargedFromProjectile,
+                                 aProjectileA-particlesFromProjectile);
     initialState2.SetNumberOfHoles(static_cast<G4int>((aProjectileA-particlesFromProjectile)/2.0));
     initialState2.SetNumberOfParticles(static_cast<G4int>((aProjectileZ-chargedFromProjectile)/2.0));

trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/src/G4PreCompoundAlpha.cc

@@ -24 +24 @@
 // ********************************************************************
 //
-//
-// $Id: G4PreCompoundAlpha.cc,v 1.6 2010/04/09 14:06:17 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-04-beta-01 $
+// $Id: G4PreCompoundAlpha.cc,v 1.7 2010/08/28 15:16:55 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03-ref-09 $
 //
 // -------------------------------------------------------------------
@@ -39 +38 @@
 // Modified:  
 // 21.08.2008 J. M. Quesada add choice of options  
+// 20.08.2010 V.Ivanchenko added G4Pow and G4PreCompoundParameters pointers
+//                         use int Z and A and cleanup
 //
 
 #include "G4PreCompoundAlpha.hh"
-
-G4ReactionProduct * G4PreCompoundAlpha::GetReactionProduct() const
-{
-  G4ReactionProduct * theReactionProduct =
-    new G4ReactionProduct(G4Alpha::AlphaDefinition());
-  theReactionProduct->SetMomentum(GetMomentum().vect());
-  theReactionProduct->SetTotalEnergy(GetMomentum().e());
-#ifdef PRECOMPOUND_TEST
-  theReactionProduct->SetCreatorModel("G4PrecompoundModel");
-#endif
-  return theReactionProduct;
-}   
-
-G4double G4PreCompoundAlpha::FactorialFactor(const G4double N, const G4double P)
-{
-  return 
-      (N-4.0)*(P-3.0)*(
-                       (((N-3.0)*(P-2.0))/2.0) *(
-                                                 (((N-2.0)*(P-1.0))/3.0) *(
-                                                                           (((N-1.0)*P)/2.0)
-                                                                           )
-                                                 )
-                       );
-}
-  
-G4double G4PreCompoundAlpha::CoalescenceFactor(const G4double A)
-{
-  return 4096.0/(A*A*A);  
+#include "G4Alpha.hh"
+
+G4PreCompoundAlpha::G4PreCompoundAlpha()
+  : G4PreCompoundIon(G4Alpha::Alpha(), &theAlphaCoulombBarrier)
+{}
+
+G4PreCompoundAlpha::~G4PreCompoundAlpha()
+{}
+
+G4double G4PreCompoundAlpha::FactorialFactor(G4int N, G4int P)
+{
+  return G4double((N-4)*(P-3)*(N-3)*(P-2)*(N-2)*(P-1)*(N-1)*P)/12.0;
+}
+  
+G4double G4PreCompoundAlpha::CoalescenceFactor(G4int A)
+{
+  return 4096.0/G4double(A*A*A);  
 }    
 
-G4double G4PreCompoundAlpha::GetRj(const G4int NumberParticles, const G4int NumberCharged)
+G4double G4PreCompoundAlpha::GetRj(G4int nParticles, G4int nCharged)
 {
   G4double rj = 0.0;
-  G4double denominator = NumberParticles*(NumberParticles-1)*(NumberParticles-2)*(NumberParticles-3);
-  if(NumberCharged >=2 && (NumberParticles-NumberCharged) >=2 ) {
-    rj = 6.0*static_cast<G4double>(NumberCharged*(NumberCharged-1)*(NumberParticles-NumberCharged)*
-                                   (NumberParticles-NumberCharged-1))/static_cast<G4double>(denominator);  
+  if(nCharged >=2 && (nParticles-nCharged) >=2 ) {
+    G4double denominator = 
+      G4double(nParticles*(nParticles-1)*(nParticles-2)*(nParticles-3));
+    rj = 6.0*nCharged*(nCharged-1)*(nParticles-nCharged)*(nParticles-nCharged-1)
+      /denominator;  
   }
   return rj;
 }
 
-////////////////////////////////////////////////////////////////////////////////////
+/////////////////////////////////////////////////////////////////////////////////
 //J. M. Quesada (Dec 2007-June 2008): New inverse reaction cross sections 
 //OPT=0 Dostrovski's parameterization
@@ -89 +80 @@
 //OPT=3,4 Kalbach's parameterization 
 // 
-G4double G4PreCompoundAlpha::CrossSection(const  G4double K)
-{
-
-  ResidualA=GetRestA();
-  ResidualZ=GetRestZ(); 
-  theA=GetA();
-  theZ=GetZ();
-  ResidualAthrd=std::pow(ResidualA,0.33333);
-  FragmentA=GetA()+GetRestA();
-  FragmentAthrd=std::pow(FragmentA,0.33333);
-
-
-  if (OPTxs==0) return GetOpt0( K);
-  else if( OPTxs==1 || OPTxs==2) return GetOpt12( K);
-  else if (OPTxs==3 || OPTxs==4)  return GetOpt34( K);
+G4double G4PreCompoundAlpha::CrossSection(G4double K)
+{
+  ResidualA = GetRestA();
+  ResidualZ = GetRestZ(); 
+  theA = GetA();
+  theZ = GetZ();
+  ResidualAthrd = ResidualA13();
+  FragmentA = theA + ResidualA;
+  FragmentAthrd = g4pow->Z13(FragmentA);
+
+  if (OPTxs==0) { return GetOpt0( K); }
+  else if( OPTxs==1 || OPTxs==2) { return GetOpt12( K); }
+  else if (OPTxs==3 || OPTxs==4) { return GetOpt34( K); }
   else{
     std::ostringstream errOs;
@@ -112 +101 @@
 }
 
-// *********************** OPT=0 : Dostrovski's cross section  *****************************
-
-G4double G4PreCompoundAlpha::GetOpt0(const  G4double K)
-{
-  const G4double r0 = G4PreCompoundParameters::GetAddress()->Getr0();
-  // cross section is now given in mb (r0 is in mm) for the sake of consistency
-  //with the rest of the options
-  return 1.e+25*pi*(r0*ResidualAthrd)*(r0*ResidualAthrd)*GetAlpha()*(1.+GetBeta()/K);
-}
-//
-//----------------
-//
 G4double G4PreCompoundAlpha::GetAlpha()
 {
   G4double C = 0.0;
-  G4double aZ = GetZ() + GetRestZ();
+  G4int aZ = theZ + ResidualZ;
   if (aZ <= 30) 
     {
@@ -134 +111 @@
   else if (aZ <= 50) 
     {
-      C = 0.1 + -((aZ-50.)/20.)*0.02;
+      C = 0.1 - (aZ-30)*0.001;
     } 
   else if (aZ < 70) 
     {
-      C = 0.08 + -((aZ-70.)/20.)*0.02;
+      C = 0.08 - (aZ-50)*0.001;
     }
   else 
@@ -146 +123 @@
   return 1.0+C;
 }
-//
-//--------------------
-// 
-G4double G4PreCompoundAlpha::GetBeta() 
-{
-  return -GetCoulombBarrier();
-}
-//
-//********************* OPT=1,2 : Chatterjee's cross section ************************ 
+
+//
+//********************* OPT=1,2 : Chatterjee's cross section ********************
 //(fitting to cross section from Bechetti & Greenles OM potential)
 
-G4double G4PreCompoundAlpha::GetOpt12(const  G4double K)
-{
-
+G4double G4PreCompoundAlpha::GetOpt12(G4double K)
+{
   G4double Kc=K;
 
   // JMQ xsec is set constant above limit of validity
-  if (K>50) Kc=50;
+  if (K > 50*MeV) { Kc = 50*MeV; }
 
   G4double landa ,mu ,nu ,p , Ec,q,r,ji,xs;
@@ -182 +152 @@
   p = p0 + p1/Ec + p2/(Ec*Ec);
   landa = landa0*ResidualA + landa1;
-  mu = mu0*std::pow(ResidualA,mu1);
-  nu = std::pow(ResidualA,mu1)*(nu0 + nu1*Ec + nu2*(Ec*Ec));
+  G4double resmu1 = g4pow->powZ(ResidualA,mu1); 
+  mu = mu0*resmu1;
+  nu = resmu1*(nu0 + nu1*Ec + nu2*(Ec*Ec));
   q = landa - nu/(Ec*Ec) - 2*p*Ec;
   r = mu + 2*nu/Ec + p*(Ec*Ec);
@@ -194 +165 @@
               
   return xs;
-
 }
 
 // *********** OPT=3,4 : Kalbach's cross sections (from PRECO code)*************
-G4double G4PreCompoundAlpha::GetOpt34(const  G4double K)
+G4double G4PreCompoundAlpha::GetOpt34(G4double K)
 // c     ** alpha from huizenga and igo
 {
-
   G4double landa, mu, nu, p , signor(1.),sig;
   G4double ec,ecsq,xnulam,etest(0.),a; 
@@ -209 +178 @@
   G4double     spill= 1.e+18;
 
-  G4double       p0 = 10.95;
+  G4double     p0 = 10.95;
   G4double     p1 = -85.2;
   G4double     p2 = 1146.;
@@ -228 +197 @@
   p = p0 + p1/ec + p2/ecsq;
   landa = landa0*ResidualA + landa1;
-  a = std::pow(ResidualA,mu1);
+  a = g4pow->powZ(ResidualA,mu1);
   mu = mu0 * a;
   nu = a* (nu0+nu1*ec+nu2*ecsq);  
   xnulam = nu / landa;
-  if (xnulam > spill) xnulam=0.;
-  if (xnulam >= flow) etest = 1.2 *std::sqrt(xnulam);
+  if (xnulam > spill) { xnulam=0.; }
+  if (xnulam >= flow) { etest = 1.2 *std::sqrt(xnulam); }
 
   a = -2.*p*ec + landa - nu/ecsq;
@@ -239 +208 @@
   ecut = 0.;
   cut = a*a - 4.*p*b;
-  if (cut > 0.) ecut = std::sqrt(cut);
+  if (cut > 0.) { ecut = std::sqrt(cut); }
   ecut = (ecut-a) / (p+p);
   ecut2 = ecut;
-//JMQ 290310 for avoiding unphysical increase below minimum (at ecut)
-//ecut<0 means that there is no cut with energy axis, i.e. xs is set to 0 bellow minimum
-//  if (cut < 0.) ecut2 = ecut - 2.;
-  if (cut < 0.) ecut2 = ecut;
-  elab = K * FragmentA / ResidualA;
+  //JMQ 290310 for avoiding unphysical increase below minimum (at ecut)
+  // ecut<0 means that there is no cut with energy axis, i.e. xs is set 
+  // to 0 bellow minimum
+  //  if (cut < 0.) ecut2 = ecut - 2.;
+  if (cut < 0.) { ecut2 = ecut; }
+  elab = K * FragmentA / G4double(ResidualA);
   sig = 0.;
   
   if (elab <= ec) { //start for E<Ec
-    if (elab > ecut2)  sig = (p*elab*elab+a*elab+b) * signor;
+    if (elab > ecut2) { sig = (p*elab*elab+a*elab+b) * signor; }
   }           //end for E<Ec
   else {           //start for E>Ec
     sig = (landa*elab+mu+nu/elab) * signor;
     geom = 0.;
-    if (xnulam < flow || elab < etest) return sig;
+    if (xnulam < flow || elab < etest) { return sig; }
     geom = std::sqrt(theA*K);
     geom = 1.23*ResidualAthrd + ra + 4.573/geom;
@@ -262 +232 @@
   }           //end for E>Ec
   return sig;
-  
-}
-
-//   ************************** end of cross sections ******************************* 
+}

trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/src/G4PreCompoundDeuteron.cc

@@ -24 +24 @@
 // ********************************************************************
 //
-// $Id: G4PreCompoundDeuteron.cc,v 1.6 2010/04/09 14:06:17 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-04-beta-01 $
+// $Id: G4PreCompoundDeuteron.cc,v 1.7 2010/08/28 15:16:55 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03-ref-09 $
 //
 // -------------------------------------------------------------------
@@ -38 +38 @@
 // Modified:  
 // 21.08.2008 J. M. Quesada add choice of options  
+// 20.08.2010 V.Ivanchenko added G4Pow and G4PreCompoundParameters pointers
+//                         use int Z and A and cleanup
 //
 
 #include "G4PreCompoundDeuteron.hh"
-
-G4ReactionProduct * G4PreCompoundDeuteron::GetReactionProduct() const
-{
-  G4ReactionProduct * theReactionProduct =
-    new G4ReactionProduct(G4Deuteron::DeuteronDefinition());
-  theReactionProduct->SetMomentum(GetMomentum().vect());
-  theReactionProduct->SetTotalEnergy(GetMomentum().e());
-#ifdef PRECOMPOUND_TEST
-  theReactionProduct->SetCreatorModel("G4PrecompoundModel");
-#endif
-  return theReactionProduct;
-}   
-
+#include "G4Deuteron.hh"
+
+G4PreCompoundDeuteron::G4PreCompoundDeuteron()
+  : G4PreCompoundIon(G4Deuteron::Deuteron(), &theDeuteronCoulombBarrier)
+{}
+
+G4PreCompoundDeuteron::~G4PreCompoundDeuteron()
+{}
  
-G4double G4PreCompoundDeuteron::FactorialFactor(const G4double N, const G4double P)
-{
-  return (N-1.0)*(N-2.0)*(P-1.0)*P/2.0;
+G4double G4PreCompoundDeuteron::FactorialFactor(G4int N, G4int P)
+{
+  return G4double((N-1)*(N-2)*(P-1)*P)/2.0;
 }
   
-G4double G4PreCompoundDeuteron::CoalescenceFactor(const G4double A)
-{
-  return 16.0/A;
+G4double G4PreCompoundDeuteron::CoalescenceFactor(G4int A)
+{
+  return 16.0/G4double(A);
 }    
 
-G4double G4PreCompoundDeuteron::GetRj(const G4int NumberParticles, const G4int NumberCharged)
+G4double G4PreCompoundDeuteron::GetRj(G4int nParticles, const G4int nCharged)
 {
   G4double rj = 0.0;
-  G4double denominator = NumberParticles*(NumberParticles-1);
-  if(NumberCharged >=1 && (NumberParticles-NumberCharged) >=1) {
-    rj = 2.0*static_cast<G4double>(NumberCharged*(NumberParticles-NumberCharged))
-      / static_cast<G4double>(denominator); 
+  if(nCharged >=1 && (nParticles-nCharged) >=1) {
+    G4double denominator = G4double(nParticles*(nParticles-1));
+    rj = 2*nCharged*(nParticles-nCharged)/denominator; 
   }
   return rj;
 }
 
-////////////////////////////////////////////////////////////////////////////////////
+////////////////////////////////////////////////////////////////////////////////
 //J. M. Quesada (Dec 2007-June 2008): New inverse reaction cross sections 
 //OPT=0 Dostrovski's parameterization
@@ -82 +78 @@
 //OPT=3,4 Kalbach's parameterization 
 // 
-G4double G4PreCompoundDeuteron::CrossSection(const  G4double K)
-{
-
-  ResidualA=GetRestA();
-  ResidualZ=GetRestZ(); 
-  theA=GetA();
-  theZ=GetZ();
-  ResidualAthrd=std::pow(ResidualA,0.33333);
-  FragmentA=GetA()+GetRestA();
-  FragmentAthrd=std::pow(FragmentA,0.33333);
-
-
-  if (OPTxs==0) return GetOpt0( K);
-  else if( OPTxs==1 || OPTxs==2) return GetOpt12( K);
-  else if (OPTxs==3 || OPTxs==4)  return GetOpt34( K);
+G4double G4PreCompoundDeuteron::CrossSection(G4double K)
+{
+  ResidualA = GetRestA();
+  ResidualZ = GetRestZ(); 
+  theA = GetA();
+  theZ = GetZ();
+  ResidualAthrd = ResidualA13();
+  FragmentA = theA + ResidualA;
+  FragmentAthrd = g4pow->Z13(FragmentA);
+
+  if (OPTxs==0) { return GetOpt0( K); }
+  else if( OPTxs==1 || OPTxs==2) { return GetOpt12( K); }
+  else if (OPTxs==3 || OPTxs==4) { return GetOpt34( K); }
   else{
     std::ostringstream errOs;
@@ -105 +99 @@
 }
 
-// *********************** OPT=0 : Dostrovski's cross section  *****************************
-
-G4double G4PreCompoundDeuteron::GetOpt0(const  G4double K)
-{
-  const G4double r0 = G4PreCompoundParameters::GetAddress()->Getr0();
-  // cross section is now given in mb (r0 is in mm) for the sake of consistency
-  //with the rest of the options
-  return 1.e+25*pi*(r0*ResidualAthrd)*(r0*ResidualAthrd)*GetAlpha()*(1.+GetBeta()/K);
-}
-//
-//---------
-//
 G4double G4PreCompoundDeuteron::GetAlpha()
 {
   G4double C = 0.0;
-  G4double aZ = GetZ() + GetRestZ();
+  G4int aZ = theZ + ResidualZ;
   if (aZ >= 70) 
     {
@@ -132 +114 @@
 }
 //
-//---------
-//
-G4double G4PreCompoundDeuteron::GetBeta() 
-{
-  return -GetCoulombBarrier();
-}
-//
-//********************* OPT=1,2 : Chatterjee's cross section ************************ 
+//********************* OPT=1,2 : Chatterjee's cross section ********************
 //(fitting to cross section from Bechetti & Greenles OM potential)
 
-G4double G4PreCompoundDeuteron::GetOpt12(const  G4double K)
-{
-
-  G4double Kc=K;
-
-// JMQ xsec is set constat above limit of validity
-  if (K>50) Kc=50;
+G4double G4PreCompoundDeuteron::GetOpt12(G4double K)
+{
+  G4double Kc = K;
+
+  // JMQ xsec is set constat above limit of validity
+  if (K > 50*MeV) { Kc = 50*MeV; }
 
   G4double landa ,mu ,nu ,p , Ec,q,r,ji,xs;
-  //G4double Eo(0),epsilon1(0),epsilon2(0),discri(0);
-
  
   G4double    p0 = -38.21;
@@ -165 +137 @@
   G4double    nu2 = -5.924;  
   G4double    delta=1.2;            
- 
 
   Ec = 1.44*theZ*ResidualZ/(1.5*ResidualAthrd+delta);
   p = p0 + p1/Ec + p2/(Ec*Ec);
   landa = landa0*ResidualA + landa1;
-  mu = mu0*std::pow(ResidualA,mu1);
-  nu = std::pow(ResidualA,mu1)*(nu0 + nu1*Ec + nu2*(Ec*Ec));
+  G4double resmu1 = g4pow->powZ(ResidualA,mu1); 
+  mu = mu0*resmu1;
+  nu = resmu1*(nu0 + nu1*Ec + nu2*(Ec*Ec));
   q = landa - nu/(Ec*Ec) - 2*p*Ec;
   r = mu + 2*nu/Ec + p*(Ec*Ec);
@@ -182 +154 @@
               
   return xs;
-
 }
 
 // *********** OPT=3,4 : Kalbach's cross sections (from PRECO code)*************
-G4double G4PreCompoundDeuteron::GetOpt34(const  G4double K)
+G4double G4PreCompoundDeuteron::GetOpt34(G4double K)
 //     ** d from o.m. of perey and perey
 {
@@ -196 +167 @@
   G4double     flow = 1.e-18;
   G4double     spill= 1.e+18;
-
- 
 
   G4double     p0 = 0.798;
@@ -210 +179 @@
   G4double     nu2 = -3.592;      
  
-
   G4double     ra=0.80;
         
@@ -219 +187 @@
   p = p0 + p1/ec + p2/ecsq;
   landa = landa0*ResidualA + landa1;
-  a = std::pow(ResidualA,mu1);
+  a = g4pow->powZ(ResidualA,mu1);
   mu = mu0 * a;
   nu = a* (nu0+nu1*ec+nu2*ecsq);  
   xnulam = nu / landa;
-  if (xnulam > spill) xnulam=0.;
-  if (xnulam >= flow) etest = 1.2 *std::sqrt(xnulam);
+  if (xnulam > spill) { xnulam=0.; }
+  if (xnulam >= flow) { etest = 1.2 *std::sqrt(xnulam); }
 
   a = -2.*p*ec + landa - nu/ecsq;
@@ -230 +198 @@
   ecut = 0.;
   cut = a*a - 4.*p*b;
-  if (cut > 0.) ecut = std::sqrt(cut);
+  if (cut > 0.) { ecut = std::sqrt(cut); }
   ecut = (ecut-a) / (p+p);
   ecut2 = ecut;
-//JMQ 290310 for avoiding unphysical increase below minimum (at ecut)
-//ecut<0 means that there is no cut with energy axis, i.e. xs is set to 0 bellow minimum
-//  if (cut < 0.) ecut2 = ecut - 2.;
-  if (cut < 0.) ecut2 = ecut;
-  elab = K * FragmentA / ResidualA;
+  //JMQ 290310 for avoiding unphysical increase below minimum (at ecut)
+  //ecut<0 means that there is no cut with energy axis, i.e. xs is set 
+  //to 0 bellow minimum
+  //  if (cut < 0.) ecut2 = ecut - 2.;
+  if (cut < 0.) { ecut2 = ecut; }
+  elab = K * FragmentA / G4double(ResidualA);
   sig = 0.;
 
   if (elab <= ec) { //start for E<Ec
-    if (elab > ecut2)  sig = (p*elab*elab+a*elab+b) * signor;
+    if (elab > ecut2) { sig = (p*elab*elab+a*elab+b) * signor; }
   }           //end for E<Ec
   else {           //start for E>Ec
     sig = (landa*elab+mu+nu/elab) * signor;
     geom = 0.;
-    if (xnulam < flow || elab < etest) return sig;
+    if (xnulam < flow || elab < etest) { return sig; }
     geom = std::sqrt(theA*K);
     geom = 1.23*ResidualAthrd + ra + 4.573/geom;
@@ -253 +222 @@
   }           //end for E>Ec
   return sig;
-
-}
-
-//   ************************** end of cross sections ******************************* 
-
-
-
-
+}

trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/src/G4PreCompoundEmission.cc

@@ -24 +24 @@
 // ********************************************************************
 //
-//
-// $Id: G4PreCompoundEmission.cc,v 1.28 2010/02/25 10:27:36 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-04-beta-01 $
+// $Id: G4PreCompoundEmission.cc,v 1.32 2010/09/01 15:11:10 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03-ref-09 $
 //
 // -------------------------------------------------------------------
@@ -41 +40 @@
 // 19.01.2010 V.Ivanchenko simplified computation of parameter an, sample cosTheta 
 //                         instead of theta; protect all calls to sqrt 
+// 20.08.2010 V.Ivanchenko added G4Pow and G4PreCompoundParameters pointers
+//                         use int Z and A and cleanup
 //
 
 #include "G4PreCompoundEmission.hh"
 #include "G4PreCompoundParameters.hh"
-
 #include "G4PreCompoundEmissionFactory.hh"
 #include "G4HETCEmissionFactory.hh"
-
-const G4PreCompoundEmission & 
-G4PreCompoundEmission::operator=(const G4PreCompoundEmission &)
-{
-  throw G4HadronicException(__FILE__, __LINE__,
-                            "G4PreCompoundEmission::operator= meant to not be accessable");
-  return *this;
-}
-
-
-G4bool G4PreCompoundEmission::operator==(const G4PreCompoundEmission &) const
-{
-  return false;
-}
-
-G4bool G4PreCompoundEmission::operator!=(const G4PreCompoundEmission &) const
-{
-  return true;
-}
+#include "G4HadronicException.hh"
+#include "G4Pow.hh"
+#include "Randomize.hh"
 
 G4PreCompoundEmission::G4PreCompoundEmission()
 {
   theFragmentsFactory = new G4PreCompoundEmissionFactory();
-  theFragmentsVector = new G4PreCompoundFragmentVector(theFragmentsFactory->GetFragmentVector());
+  theFragmentsVector = 
+    new G4PreCompoundFragmentVector(theFragmentsFactory->GetFragmentVector());
+  g4pow = G4Pow::GetInstance();
+  theParameters = G4PreCompoundParameters::GetAddress();
 }
 
 G4PreCompoundEmission::~G4PreCompoundEmission()
 {
-  if (theFragmentsFactory) delete theFragmentsFactory;
-  if (theFragmentsVector) delete theFragmentsVector;
+  if (theFragmentsFactory) { delete theFragmentsFactory; }
+  if (theFragmentsVector)  { delete theFragmentsVector; }
 }
 
 void G4PreCompoundEmission::SetDefaultModel()
 {
-  if (theFragmentsFactory) delete theFragmentsFactory;
+  if (theFragmentsFactory) { delete theFragmentsFactory; }
   theFragmentsFactory = new G4PreCompoundEmissionFactory();
   if (theFragmentsVector) 
@@ -90 +77 @@
   else 
     {
-      theFragmentsVector = new G4PreCompoundFragmentVector(theFragmentsFactory->GetFragmentVector());
-    }
-  theFragmentsVector->ResetStage();
+      theFragmentsVector = 
+        new G4PreCompoundFragmentVector(theFragmentsFactory->GetFragmentVector());
+    }
   return;
 }
@@ -106 +93 @@
   else 
     {
-      theFragmentsVector = new G4PreCompoundFragmentVector(theFragmentsFactory->GetFragmentVector());
-    }
-  theFragmentsVector->ResetStage();
+      theFragmentsVector = 
+        new G4PreCompoundFragmentVector(theFragmentsFactory->GetFragmentVector());
+    }
   return;
 }
@@ -114 +101 @@
 G4ReactionProduct * G4PreCompoundEmission::PerformEmission(G4Fragment & aFragment)
 {
-#ifdef debug
-  G4Fragment InitialState(aFragment);
-#endif
   // Choose a Fragment for emission
-  G4VPreCompoundFragment * theFragment = theFragmentsVector->ChooseFragment();
-  if (theFragment == 0)
-    {
-      G4cerr <<  "G4PreCompoundEmission::PerformEmission : I couldn't choose a fragment\n"
+  G4VPreCompoundFragment * thePreFragment = theFragmentsVector->ChooseFragment();
+  if (thePreFragment == 0)
+    {
+      G4cout <<  "G4PreCompoundEmission::PerformEmission : I couldn't choose a fragment\n"
              << "while trying to de-excite\n" 
-             << aFragment << '\n';
+             << aFragment << G4endl;
       throw G4HadronicException(__FILE__, __LINE__, "");
     }
   // Kinetic Energy of emitted fragment
-  G4double KineticEnergyOfEmittedFragment = theFragment->GetKineticEnergy(aFragment);
+  G4double kinEnergyOfEmittedFragment = thePreFragment->GetKineticEnergy(aFragment);
+  if(kinEnergyOfEmittedFragment < 0.0) { kinEnergyOfEmittedFragment = 0.0; }
   
   // Calculate the fragment momentum (three vector)
-  G4ThreeVector momentum = AngularDistribution(theFragment,aFragment,KineticEnergyOfEmittedFragment);
+  AngularDistribution(thePreFragment,aFragment,kinEnergyOfEmittedFragment);
   
   // Mass of emittef fragment
-  G4double EmittedMass = theFragment->GetNuclearMass();
+  G4double EmittedMass = thePreFragment->GetNuclearMass();
   
   // Now we can calculate the four momentum 
   // both options are valid and give the same result but 2nd one is faster
-  // G4LorentzVector EmittedMomentum(momentum,std::sqrt(momentum.mag2()+EmittedMass*EmittedMass));
-  G4LorentzVector EmittedMomentum(momentum,EmittedMass+KineticEnergyOfEmittedFragment);
+  G4LorentzVector Emitted4Momentum(theFinalMomentum,
+                                   EmittedMass + kinEnergyOfEmittedFragment);
     
   // Perform Lorentz boost
-  EmittedMomentum.boost(aFragment.GetMomentum().boostVector());  
+  G4LorentzVector Rest4Momentum = aFragment.GetMomentum();
+  Emitted4Momentum.boost(Rest4Momentum.boostVector());  
 
   // Set emitted fragment momentum
-  theFragment->SetMomentum(EmittedMomentum);    
-
+  thePreFragment->SetMomentum(Emitted4Momentum);        
 
   // NOW THE RESIDUAL NUCLEUS
   // ------------------------
-    
-  // Now the residual nucleus. 
-  // The energy conservation says that 
-  G4double ResidualEcm = 
-    //    aFragment.GetGroundStateMass() + aFragment.GetExcitationEnergy() // initial energy in cm
-    aFragment.GetMomentum().m()
-    - (EmittedMass+KineticEnergyOfEmittedFragment); 
-
-  // Then the four momentum for residual is 
-  G4LorentzVector RestMomentum(-momentum,ResidualEcm);
-  // This could save a Lorentz boost
-  // G4LorentzVector RestMomentum2(aFragment.GetMomentum()-EmittedMomentum);
-
-  // Just for test
-  // Excitation energy
-  //  G4double anU = ResidualEcm - theFragment->GetRestNuclearMass();
-  // This is equivalent
-  //  G4double anU = theFragment->GetMaximalKineticEnergy() - KineticEnergyOfEmittedFragment + 
-  //    theFragment->GetCoulombBarrier();
-    
-  // check that Excitation energy is >= 0
-  G4double anU = RestMomentum.m()-theFragment->GetRestNuclearMass();
-  if (anU < 0.0) {
-    throw G4HadronicException(__FILE__, __LINE__, 
-                              "G4PreCompoundModel::DeExcite: Excitation energy less than 0!");
-  }
-     
+
+  Rest4Momentum -= Emitted4Momentum;
+    
   // Update nucleus parameters:
   // --------------------------
@@ -181 +142 @@
   // Number of excitons
   aFragment.SetNumberOfParticles(aFragment.GetNumberOfParticles()-
-                                 static_cast<G4int>(theFragment->GetA()));
+                                 thePreFragment->GetA());
   // Number of charges
   aFragment.SetNumberOfCharged(aFragment.GetNumberOfCharged()-
-                               static_cast<G4int>(theFragment->GetZ()));
-    
-  // Atomic number
-  aFragment.SetA(theFragment->GetRestA());
-    
-  // Charge
-  aFragment.SetZ(theFragment->GetRestZ());
-
-    
-  // Perform Lorentz boosts
-  RestMomentum.boost(aFragment.GetMomentum().boostVector());
-
-  // Update nucleus momentum
-  aFragment.SetMomentum(RestMomentum);
+                               thePreFragment->GetZ());
+    
+  // Z and A
+  aFragment.SetZandA_asInt(thePreFragment->GetRestZ(),
+                           thePreFragment->GetRestA());
+    
+  // Update nucleus momentum 
+  // A check on consistence of Z, A, and mass will be performed
+  aFragment.SetMomentum(Rest4Momentum);
         
   // Create a G4ReactionProduct 
-  G4ReactionProduct * MyRP = theFragment->GetReactionProduct();
-#ifdef PRECOMPOUND_TEST
-  MyRP->SetCreatorModel("G4PreCompoundModel");
-#endif
-#ifdef debug
-  CheckConservation(InitialState,aFragment,MyRP);
-#endif
+  G4ReactionProduct * MyRP = thePreFragment->GetReactionProduct();
+
+  //G4cout << "G4PreCompoundEmission::Fragment emitted" << G4endl;
+  //G4cout << thePreFragment << G4endl;
+
   return MyRP;
 }
 
-G4ThreeVector 
-G4PreCompoundEmission::AngularDistribution(G4VPreCompoundFragment * theFragment,
+void
+G4PreCompoundEmission::AngularDistribution(G4VPreCompoundFragment* thePreFragment,
                                            const G4Fragment& aFragment,
-                                           const G4double kinEnergyOfEmittedFrag) const
-{
-  G4double p = aFragment.GetNumberOfParticles();
-  G4double h = aFragment.GetNumberOfHoles();
+                                           G4double ekin) 
+{
+  G4int p = aFragment.GetNumberOfParticles();
+  G4int h = aFragment.GetNumberOfHoles();
   G4double U = aFragment.GetExcitationEnergy();
 
-  G4double ekin = std::max(0.0, kinEnergyOfEmittedFrag);
-        
   // Emission particle separation energy
-  G4double Bemission = theFragment->GetBindingEnergy();
+  G4double Bemission = thePreFragment->GetBindingEnergy();
 
   // Fermi energy
-  G4double Ef = G4PreCompoundParameters::GetAddress()->GetFermiEnergy();
+  G4double Ef = theParameters->GetFermiEnergy();
         
   //
@@ -231 +183 @@
   //  G4double g = (6.0/pi2)*aFragment.GetA()*
 
-  G4double g = (6.0/pi2)*aFragment.GetA()
-    *G4PreCompoundParameters::GetAddress()->GetLevelDensity();
+  G4double g = (6.0/pi2)*aFragment.GetA_asInt()*theParameters->GetLevelDensity();
         
   // Average exciton energy relative to bottom of nuclear well
-  G4double Eav = 2.0*p*(p+1.0)/((p+h)*g);
+  G4double Eav = 2*p*(p+1)/((p+h)*g);
         
   // Excitation energy relative to the Fermi Level
@@ -241 +192 @@
   //  G4double Uf = U - KineticEnergyOfEmittedFragment - Bemission;
 
-  G4double w_num = rho(p+1,h,g,Uf,Ef);
-  G4double w_den = rho(p,h,g,Uf,Ef);
+  G4double w_num = rho(p+1, h, g, Uf, Ef);
+  G4double w_den = rho(p,   h, g, Uf, Ef);
   if (w_num > 0.0 && w_den > 0.0)
     {
@@ -253 +204 @@
     }
   
-
   // VI + JMQ 19/01/2010 update computation of the parameter an
   //
@@ -262 +212 @@
     G4double zeta = std::max(1.0,9.3/std::sqrt(ekin/MeV));
   
-    an = 3.0*std::sqrt((ProjEnergy+Ef)*Eeff)/(zeta*Eav);
+    // This should be the projectile energy. If I would know which is 
+    // the projectile (proton, neutron) I could remove the binding energy. 
+    // But, what happens if INC precedes precompound? This approximation
+    // seems to work well enough
+    G4double ProjEnergy = aFragment.GetExcitationEnergy();
+
+    an = 3*std::sqrt((ProjEnergy+Ef)*Eeff)/(zeta*Eav);
 
     G4int ne = aFragment.GetNumberOfExcitons() - 1;
@@ -283 +239 @@
   }  
 
-  G4double phi = twopi*G4UniformRand();
+  G4double phi = CLHEP::twopi*G4UniformRand();
   
   // Calculate the momentum magnitude of emitted fragment       
-  G4double pmag = std::sqrt(ekin*(ekin + 2.0*theFragment->GetNuclearMass()));
+  G4double pmag = std::sqrt(ekin*(ekin + 2.0*thePreFragment->GetNuclearMass()));
   
   G4double sint = std::sqrt((1.0-cost)*(1.0+cost));
 
-  G4ThreeVector momentum(pmag*std::cos(phi)*sint,pmag*std::sin(phi)*sint,pmag*cost);
+  theFinalMomentum.set(pmag*std::cos(phi)*sint,pmag*std::sin(phi)*sint,pmag*cost);
+
   // theta is the angle wrt the incident direction
-  momentum.rotateUz(theIncidentDirection);
-
-  return momentum;
-}
-
-G4double G4PreCompoundEmission::rho(const G4double p, const G4double h, const G4double g, 
-                                    const G4double E, const G4double Ef) const
+  G4ThreeVector theIncidentDirection = aFragment.GetMomentum().vect().unit();
+  theFinalMomentum.rotateUz(theIncidentDirection);
+}
+
+G4double G4PreCompoundEmission::rho(G4int p, G4int h, G4double g, 
+                                    G4double E, G4double Ef) const
 {       
   // 25.02.2010 V.Ivanchenko added more protections
@@ -306 +262 @@
   if ( E - Aph < 0.0) { return 0.0; }
   
-  G4double logConst =  (p+h)*std::log(g) - logfactorial(p+h-1) - logfactorial(p) - logfactorial(h);
+  G4double logConst =  (p+h)*std::log(g) 
+    - g4pow->logfactorial(p+h-1) - g4pow->logfactorial(p) - g4pow->logfactorial(h);
 
   // initialise values using j=0
@@ -312 +269 @@
   G4double t1=1;
   G4double t2=1;
-  G4double logt3=(p+h-1) * std::log(E-Aph) + logConst;
+  G4double logt3 = (p+h-1) * std::log(E-Aph) + logConst;
   const G4double logmax = 200.;
   if(logt3 > logmax) { logt3 = logmax; }
@@ -333 +290 @@
   return tot;
 }
-
-G4double G4PreCompoundEmission::factorial(G4double a) const
-{
-  // Values of factorial function from 0 to 60
-  const G4int factablesize = 61;
-  static const G4double fact[factablesize] = 
-    {
-      1.0, // 0!
-      1.0, // 1!
-      2.0, // 2!
-      6.0, // 3!
-      24.0, // 4!
-      120.0, // 5!
-      720.0, // 6!
-      5040.0, // 7!
-      40320.0, // 8!
-      362880.0, // 9!
-      3628800.0, // 10!
-      39916800.0, // 11!
-      479001600.0, // 12!
-      6227020800.0, // 13!
-      87178291200.0, // 14!
-      1307674368000.0, // 15!
-      20922789888000.0, // 16!
-      355687428096000.0, // 17!
-      6402373705728000.0, // 18!
-      121645100408832000.0, // 19!
-      2432902008176640000.0, // 20!
-      51090942171709440000.0, // 21!
-      1124000727777607680000.0, // 22!
-      25852016738884976640000.0, // 23!
-      620448401733239439360000.0, // 24!
-      15511210043330985984000000.0, // 25!
-      403291461126605635584000000.0, // 26!
-      10888869450418352160768000000.0, // 27!
-      304888344611713860501504000000.0, // 28!
-      8841761993739701954543616000000.0, // 29!
-      265252859812191058636308480000000.0, // 30!
-      8222838654177922817725562880000000.0, // 31!
-      263130836933693530167218012160000000.0, // 32!
-      8683317618811886495518194401280000000.0, // 33!
-      295232799039604140847618609643520000000.0, // 34!
-      10333147966386144929666651337523200000000.0, // 35!
-      371993326789901217467999448150835200000000.0, // 36!
-      13763753091226345046315979581580902400000000.0, // 37!
-      523022617466601111760007224100074291200000000.0, // 38!
-      20397882081197443358640281739902897356800000000.0, // 39!
-      815915283247897734345611269596115894272000000000.0, // 40!
-      33452526613163807108170062053440751665152000000000.0, // 41!
-      1405006117752879898543142606244511569936384000000000.0, // 42!
-      60415263063373835637355132068513997507264512000000000.0, // 43!
-      2658271574788448768043625811014615890319638528000000000.0, // 44!
-      119622220865480194561963161495657715064383733760000000000.0, // 45!
-      5502622159812088949850305428800254892961651752960000000000.0, // 46!
-      258623241511168180642964355153611979969197632389120000000000.0, // 47!
-      12413915592536072670862289047373375038521486354677760000000000.0, // 48!
-      608281864034267560872252163321295376887552831379210240000000000.0, // 49!
-      30414093201713378043612608166064768844377641568960512000000000000.0, // 50!
-      1551118753287382280224243016469303211063259720016986112000000000000.0, // 51!
-      80658175170943878571660636856403766975289505440883277824000000000000.0, // 52!
-      4274883284060025564298013753389399649690343788366813724672000000000000.0, // 53!
-      230843697339241380472092742683027581083278564571807941132288000000000000.0, // 54!
-      12696403353658275925965100847566516959580321051449436762275840000000000000.0, // 55!
-      710998587804863451854045647463724949736497978881168458687447040000000000000.0, // 56!
-      40526919504877216755680601905432322134980384796226602145184481280000000000000.0, // 57!
-      2350561331282878571829474910515074683828862318181142924420699914240000000000000.0, // 58!
-      138683118545689835737939019720389406345902876772687432540821294940160000000000000.0, // 59!
-      8320987112741390144276341183223364380754172606361245952449277696409600000000000000.0  // 60!
-    };
-  //    fact[0] = 1;
-  //    for (G4int n = 1; n < 21; n++) {
-  //      fact[n] = fact[n-1]*static_cast<G4double>(n); 
-  //    }
-  G4double result(0.0);
-  G4int ia = static_cast<G4int>(a);
-  if (ia < factablesize) 
-    {
-      result = fact[ia];
-    }
-  else
-    {
-      result = fact[factablesize-1];
-      for (G4int n = factablesize; n < ia+1; ++n)
-        {
-          result *= static_cast<G4double>(n);
-        }
-    }
-    
-    return result;
-}
-G4double G4PreCompoundEmission::logfactorial(G4double a) const
-{
-  // Values of logs of factorial function from 0 to 60
-
-  G4double result(0.0);
-  const G4int factablesize = 61;
-  const G4double halfLn2pi = 0.918938533;      // 0.5 log(2* pi)
-  static G4double logfact[factablesize];
-  static bool needinit=true;
-  
-  if (needinit)
-  {
-      needinit=false;
-      for ( G4int n=0; n < factablesize; ++n)
-      {
-         logfact[n]=std::log(factorial(n));
-      }
-  }
-
-  G4int ia = static_cast<G4int>(a);
-  if (ia < factablesize)
-  {
-      result = logfact[ia];
-  } else {
-      result = (ia+0.5)*std::log(G4double(ia)) - ia + halfLn2pi;
-  }
-   
-  return result;
-}
-
-#ifdef debug
-void G4PreCompoundEmission::CheckConservation(const G4Fragment & theInitialState,
-                                              const G4Fragment & theResidual,
-                                              G4ReactionProduct * theEmitted) const
-{
-  G4double ProductsEnergy = theEmitted->GetTotalEnergy() + theResidual.GetMomentum().e();
-  G4ThreeVector ProductsMomentum(theEmitted->GetMomentum()+theResidual.GetMomentum().vect());
-  G4int ProductsA = theEmitted->GetDefinition()->GetBaryonNumber() + theResidual.GetA();
-  G4int ProductsZ = theEmitted->GetDefinition()->GetPDGCharge() + theResidual.GetZ();
-
-  if (ProductsA != theInitialState.GetA()) {
-    G4cout << "!!!!!!!!!! Baryonic Number Conservation Violation !!!!!!!!!!" << G4endl;
-    G4cout << "G4PreCompoundEmission.cc: Barionic Number Conservation"
-           << G4endl; 
-    G4cout << "Initial A = " << theInitialState.GetA() 
-           << "   Fragments A = " << ProductsA << "   Diference --> " 
-           << theInitialState.GetA() - ProductsA << G4endl;
-  }
-  if (ProductsZ != theInitialState.GetZ()) {
-    G4cout << "!!!!!!!!!! Charge Conservation Violation !!!!!!!!!!" << G4endl;
-    G4cout << "G4PreCompoundEmission.cc: Charge Conservation test"
-           << G4endl; 
-    G4cout << "Initial Z = " << theInitialState.GetZ() 
-           << "   Fragments Z = " << ProductsZ << "   Diference --> " 
-           << theInitialState.GetZ() - ProductsZ << G4endl;
-  }
-  if (std::abs(ProductsEnergy-theInitialState.GetMomentum().e()) > 10.0*eV) {
-    G4cout << "!!!!!!!!!! Energy Conservation Violation !!!!!!!!!!" << G4endl;
-    G4cout << "G4PreCompoundEmission.cc: Energy Conservation test" 
-           << G4endl; 
-    G4cout << "Initial E = " << theInitialState.GetMomentum().e()/MeV << " MeV"
-           << "   Fragments E = " << ProductsEnergy/MeV  << " MeV   Diference --> " 
-           << (theInitialState.GetMomentum().e() - ProductsEnergy)/MeV << " MeV" << G4endl;
-  } 
-  if (std::abs(ProductsMomentum.x()-theInitialState.GetMomentum().x()) > 10.0*eV || 
-      std::abs(ProductsMomentum.y()-theInitialState.GetMomentum().y()) > 10.0*eV ||
-      std::abs(ProductsMomentum.z()-theInitialState.GetMomentum().z()) > 10.0*eV) {
-    G4cout << "!!!!!!!!!! Momentum Conservation Violation !!!!!!!!!!" << G4endl;
-    G4cout << "G4PreCompoundEmission.cc: Momentum Conservation test" 
-           << G4endl; 
-    G4cout << "Initial P = " << theInitialState.GetMomentum().vect() << " MeV"
-           << "   Fragments P = " << ProductsMomentum  << " MeV   Diference --> " 
-           << theInitialState.GetMomentum().vect() - ProductsMomentum << " MeV" << G4endl;
-  }
-  return;
-}
-
-#endif

trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/src/G4PreCompoundEmissionFactory.cc

@@ -24 +24 @@
 // ********************************************************************
 //
+// $Id: G4PreCompoundEmissionFactory.cc,v 1.5 2010/08/28 15:16:55 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03-ref-09 $
+//
+
 #include "G4PreCompoundEmissionFactory.hh"
 
@@ -33 +37 @@
 #include "G4PreCompoundAlpha.hh"
 
+G4PreCompoundEmissionFactory::G4PreCompoundEmissionFactory()
+{}
 
-const G4PreCompoundEmissionFactory & G4PreCompoundEmissionFactory::
-operator=(const G4PreCompoundEmissionFactory & )
-{
-  throw G4HadronicException(__FILE__, __LINE__, "G4PreCompoundEmissionFactory::operator= meant to not be accessable.");
-  return *this;
-}
-
-G4bool G4PreCompoundEmissionFactory::
-operator==(const G4PreCompoundEmissionFactory & ) const
-{
-  throw G4HadronicException(__FILE__, __LINE__, "G4PreCompoundEmissionFactory::operator== meant to not be accessable.");
-  return false;
-}
-
-G4bool G4PreCompoundEmissionFactory::
-operator!=(const G4PreCompoundEmissionFactory & ) const
-{
-  throw G4HadronicException(__FILE__, __LINE__, "G4PreCompoundEmissionFactory::operator!= meant to not be accessable.");
-  return true;
-}
-
+G4PreCompoundEmissionFactory::~G4PreCompoundEmissionFactory()
+{}
 
 std::vector<G4VPreCompoundFragment*> *  G4PreCompoundEmissionFactory::

trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/src/G4PreCompoundFragment.cc

@@ -24 +24 @@
 // ********************************************************************
 //
-// $Id: G4PreCompoundFragment.cc,v 1.8 2009/02/10 16:01:37 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-04-beta-01 $
+// $Id: G4PreCompoundFragment.cc,v 1.9 2010/08/28 15:16:55 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03-ref-09 $
 //
 // J. M. Quesada (August 2008).  
 // Based  on previous work by V. Lara
-// JMQ (06 September 2008) Also external choice has been added for:
-//                      - superimposed Coulomb barrier (if useSICB=true) 
 //
+// Modified:
+// 06.09.2008 JMQ Also external choice has been added for:
+//               - superimposed Coulomb barrier (if useSICB=true) 
+// 20.08.2010 V.Ivanchenko cleanup
+//
+
 #include "G4PreCompoundFragment.hh"
 
 G4PreCompoundFragment::
-G4PreCompoundFragment(const G4PreCompoundFragment &right) :
-  G4VPreCompoundFragment(right)
+G4PreCompoundFragment(const G4ParticleDefinition* part,
+                      G4VCoulombBarrier* aCoulombBarrier)
+  : G4VPreCompoundFragment(part,aCoulombBarrier)
 {}
 
-
-G4PreCompoundFragment::
-G4PreCompoundFragment(const G4double anA,
-                      const G4double aZ, 
-                      G4VCoulombBarrier* aCoulombBarrier,
-                      const G4String & aName):
-  G4VPreCompoundFragment(anA,aZ,aCoulombBarrier,aName)
-{
-}
-
-
-
 G4PreCompoundFragment::~G4PreCompoundFragment()
-{
-}
-
-
-const G4PreCompoundFragment & G4PreCompoundFragment::
-operator= (const G4PreCompoundFragment & right)
-{
-  if (&right != this) this->G4VPreCompoundFragment::operator=(right);
-  return *this;
-}
-
-G4int G4PreCompoundFragment::operator==(const G4PreCompoundFragment & right) const
-{
-  return G4VPreCompoundFragment::operator==(right);
-}
-
-G4int G4PreCompoundFragment::operator!=(const G4PreCompoundFragment & right) const
-{
-  return G4VPreCompoundFragment::operator!=(right);
-}
-
+{}
 
 G4double G4PreCompoundFragment::
 CalcEmissionProbability(const G4Fragment & aFragment)
 {
-// If  theCoulombBarrier effect is included in the emission probabilities
-//if (GetMaximalKineticEnergy() <= 0.0) 
-  G4double limit;
-  if(OPTxs==0 ||  useSICB) limit= theCoulombBarrier;
-  else limit=0.;
+  //G4cout << theCoulombBarrier << "  " << GetMaximalKineticEnergy() << G4endl;
+  // If  theCoulombBarrier effect is included in the emission probabilities
+  //if (GetMaximalKineticEnergy() <= 0.0) 
+  G4double limit = 0.0;
+  if(OPTxs==0 ||  useSICB) { limit = theCoulombBarrier; }
   if (GetMaximalKineticEnergy() <= limit) 
     {
       theEmissionProbability = 0.0;
       return 0.0;
-  }    
-// If  theCoulombBarrier effect is included in the emission probabilities
-//  G4double LowerLimit = 0.;
-// Coulomb barrier is the lower limit 
-// of integration over kinetic energy
+    }    
+  // If  theCoulombBarrier effect is included in the emission probabilities
+  //  G4double LowerLimit = 0.;
+  // Coulomb barrier is the lower limit 
+  // of integration over kinetic energy
   G4double LowerLimit = limit;
 
-// Excitation energy of nucleus after fragment emission is the upper limit of integration over kinetic energy
+  // Excitation energy of nucleus after fragment emission is the upper 
+  //limit of integration over kinetic energy
   G4double UpperLimit = GetMaximalKineticEnergy();
   
   theEmissionProbability = 
     IntegrateEmissionProbability(LowerLimit,UpperLimit,aFragment);
+  /*
+  G4cout << "## G4PreCompoundFragment::CalcEmisProb "
+         << "Z= " << aFragment.GetZ_asInt() 
+         << " A= " << aFragment.GetA_asInt()
+         << " Elow= " << LowerLimit/MeV
+         << " Eup= " << UpperLimit/MeV
+         << " prob= " << theEmissionProbability
+         << G4endl;
+  */
   return theEmissionProbability;
 }
 
 G4double G4PreCompoundFragment::
-IntegrateEmissionProbability(const G4double & Low, const G4double & Up,
+IntegrateEmissionProbability(G4double Low, G4double Up,
                              const G4Fragment & aFragment)
 {
@@ -134 +117 @@
   G4double Total = 0.0;
 
-
-  for (G4int i = 0; i < N; i++) 
+  for (G4int i = 0; i < N; ++i) 
     {
-      G4double KineticE = ((Up-Low)*x[i]+(Up+Low))/2.0;
+      G4double KineticE = 0.5*((Up-Low)*x[i]+(Up+Low));
       Total += w[i]*ProbabilityDistributionFunction(KineticE, aFragment);
     }
-  Total  *= (Up-Low)/2.0;
+  Total  *= 0.5*(Up-Low);
   return Total;
 }
-
-
-
 
 G4double G4PreCompoundFragment::
 GetKineticEnergy(const G4Fragment & aFragment) 
 {
+  //let's keep this way for consistency with CalcEmissionProbability method
+  G4double V = 0.0;
+  if(OPTxs==0 || useSICB) { V = theCoulombBarrier; }
 
-//      G4double V = this->GetCoulombBarrier();// alternative way for accessing the Coulomb barrier
-//                                             //should be equivalent (in fact it is)
-  G4double V;
-  if(OPTxs==0 || useSICB) V= theCoulombBarrier;//let's keep this way for consistency with CalcEmissionProbability method
-  else V=0.;
-
-  G4double Tmax =  GetMaximalKineticEnergy() ;  
+  G4double Tmax = GetMaximalKineticEnergy();  
+  if(Tmax < V) { return 0.0; }
   G4double T(0.0);
-  G4double NormalizedProbability(1.0);
+  G4double Probability(1.0);
+  G4double maxProbability = GetEmissionProbability();
   do 
     {
-      T =V+ G4UniformRand()*(Tmax-V);
-      NormalizedProbability = ProbabilityDistributionFunction(T,aFragment)/GetEmissionProbability();      
-    }   while (G4UniformRand() > NormalizedProbability);  
+      T = V + G4UniformRand()*(Tmax-V);
+      Probability = ProbabilityDistributionFunction(T,aFragment);      
+    }   while (maxProbability*G4UniformRand() > Probability);  
   return T;
 }

trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/src/G4PreCompoundFragmentVector.cc

@@ -24 +24 @@
 // ********************************************************************
 //
-// $Id: G4PreCompoundFragmentVector.cc,v 1.11 2009/02/10 16:01:37 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-04-beta-01 $
+// $Id: G4PreCompoundFragmentVector.cc,v 1.12 2010/08/28 15:16:55 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03-ref-09 $
 //
 // Hadronic Process: Nuclear Preequilibrium
 // by V. Lara 
+//
+// Modified:
+// 27.08.2010 V.Ivanchenko moved constructor and destructor to source, 
+//            simplify run time computations making inlined 
+//
 
 #include "G4PreCompoundFragmentVector.hh"
-#include "G4HadronicException.hh"
 
-const G4PreCompoundFragmentVector & 
-G4PreCompoundFragmentVector::
-operator=(const G4PreCompoundFragmentVector &)
+G4PreCompoundFragmentVector::G4PreCompoundFragmentVector(pcfvector * avector) 
+  : theChannels(0), nChannels(0) 
 {
-    throw G4HadronicException(__FILE__, __LINE__, "G4PreCompoundFragmentVector::operator= meant to not be accessable");
-    return *this;
+  SetVector(avector);
 }
 
+G4PreCompoundFragmentVector::~G4PreCompoundFragmentVector()
+{}
 
-G4bool G4PreCompoundFragmentVector::
-operator==(const G4PreCompoundFragmentVector &) const
+void G4PreCompoundFragmentVector::SetVector(pcfvector * avector)
 {
-    return false;
+  theChannels = avector;
+  if(theChannels) {
+    nChannels = theChannels->size();
+    probabilities.resize(nChannels);
+  }
 }
 
-G4bool G4PreCompoundFragmentVector::
-operator!=(const G4PreCompoundFragmentVector &) const
-{
-    return true;
+//for inverse cross section choice
+void G4PreCompoundFragmentVector::SetOPTxs(G4int opt)
+{    
+  for (G4int i=0; i< nChannels; ++i) { 
+    (*theChannels)[i]->SetOPTxs(opt);
+  }
 }
 
-
-
-G4double G4PreCompoundFragmentVector::
-CalculateProbabilities(const G4Fragment & aFragment)
-{
-  TotalEmissionProbability = 0.0;
-  pcfvector::iterator aChannel; 
-  for (aChannel=theChannels->begin(); aChannel != theChannels->end(); 
-       aChannel++) 
-    {
-      // Calculate emission probailities
-      // Compute total (integrated over kinetic energy) emission 
-      // probability of a fragment and
-      // Summing channel emission probabilities
-      TotalEmissionProbability += (*aChannel)->CalcEmissionProbability(aFragment);
-    }
-  return TotalEmissionProbability;
+//for superimposed Coulomb Barrier for inverse  cross sections 
+void G4PreCompoundFragmentVector::UseSICB(G4bool use)
+{    
+  for (G4int i=0; i< nChannels; ++i) { 
+    (*theChannels)[i]->UseSICB(use);
+  }
 }
 
-
-G4VPreCompoundFragment * G4PreCompoundFragmentVector::
-ChooseFragment(void)
-{
-  const G4int NumOfFrags = theChannels->size();
-  std::vector<G4double> running;
-  running.reserve(NumOfFrags);
-  
-  pcfvector::iterator i;
-  G4double accumulation = 0.0;
-  for (i = theChannels->begin(); i != theChannels->end(); ++i) {
-    accumulation += (*i)->GetEmissionProbability();
-
-    running.push_back(accumulation);
-  }
-        
-  // Choose an emission channel
-  G4double aChannel = G4UniformRand()*TotalEmissionProbability;
-  G4int ChosenChannel = -1;
-  std::vector<G4double>::iterator ich;
-  for (ich = running.begin(); ich != running.end(); ++ich) 
-    {
-      if (aChannel <= *ich) 
-        {
-#ifdef G4NO_ISO_VECDIST
-          std::vector<G4double>::difference_type n = 0;
-          std::distance(running.begin(),ich,n);
-          ChosenChannel = n;
-#else
-          ChosenChannel = std::distance(running.begin(),ich);
-#endif
-          break;
-        }
-    }
-  running.clear();
-  if (ChosenChannel < 0) 
-    {
-      G4cerr
-        << "G4PreCompoundFragmentVector::ChooseFragment: I can't determine a channel\n"
-        << "Probabilities: ";
-      for (i = theChannels->begin(); i != theChannels->end(); ++i) 
-        {
-          G4cout << (*i)->GetEmissionProbability() << "  ";
-        }
-      G4cout << '\n';
-      return 0;
-    }
-  else
-    {
-      for (i = theChannels->begin(); i != theChannels->end(); ++i) 
-        {
-          (*i)->IncrementStage();
-        }
-    }
-
-  return theChannels->operator[](ChosenChannel);
-}

trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/src/G4PreCompoundHe3.cc

@@ -24 +24 @@
 // ********************************************************************
 //
-//
-// $Id: G4PreCompoundHe3.cc,v 1.6 2010/04/09 14:06:17 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-04-beta-01 $
+// $Id: G4PreCompoundHe3.cc,v 1.7 2010/08/28 15:16:55 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03-ref-09 $
 //
 // -------------------------------------------------------------------
@@ -39 +38 @@
 // Modified:  
 // 21.08.2008 J. M. Quesada add choice of options  
+// 20.08.2010 V.Ivanchenko added G4Pow and G4PreCompoundParameters pointers
+//                         use int Z and A and cleanup
 //
  
 #include "G4PreCompoundHe3.hh"
-
-G4ReactionProduct * G4PreCompoundHe3::GetReactionProduct() const
-{
-  G4ReactionProduct * theReactionProduct =
-    new G4ReactionProduct(G4He3::He3Definition());
-  theReactionProduct->SetMomentum(GetMomentum().vect());
-  theReactionProduct->SetTotalEnergy(GetMomentum().e());
-#ifdef PRECOMPOUND_TEST
-  theReactionProduct->SetCreatorModel("G4PrecompoundModel");
-#endif
-  return theReactionProduct;
-}   
-
-G4double G4PreCompoundHe3::FactorialFactor(const G4double N, const G4double P)
-{
-  return 
-      (N-3.0)*(P-2.0)*(
-                       (((N-2.0)*(P-1.0))/2.0) *(
-                                                 (((N-1.0)*P)/3.0) 
-                                                 )
-                       );
-}
-  
-G4double G4PreCompoundHe3::CoalescenceFactor(const G4double A)
-{
-  return 243.0/(A*A);
+#include "G4He3.hh"
+
+G4PreCompoundHe3::G4PreCompoundHe3()
+  : G4PreCompoundIon(G4He3::He3(), &theHe3CoulombBarrier)
+{}
+
+G4PreCompoundHe3::~G4PreCompoundHe3()
+{}
+
+G4double G4PreCompoundHe3::FactorialFactor(G4int N, G4int P)
+{
+  return G4double((N-3)*(P-2)*(N-2)*(P-1)*(N-1)*P)/6.0; 
+}
+  
+G4double G4PreCompoundHe3::CoalescenceFactor(G4int A)
+{
+  return 243.0/G4double(A*A);
 }    
 
-G4double G4PreCompoundHe3::GetRj(const G4int NumberParticles, const G4int NumberCharged)
+G4double G4PreCompoundHe3::GetRj(G4int nParticles, G4int nCharged)
 {
   G4double rj = 0.0;
-  G4double denominator = NumberParticles*(NumberParticles-1)*(NumberParticles-2);
-  if(NumberCharged >=2 && (NumberParticles-NumberCharged) >= 1) {
-    rj = 3.0*static_cast<G4double>(NumberCharged*(NumberCharged-1)*(NumberParticles-NumberCharged))
-      / static_cast<G4double>(denominator);  
+  if(nCharged >=2 && (nParticles-nCharged) >= 1) {
+    G4double denominator = G4double(nParticles*(nParticles-1)*(nParticles-2));
+    rj = G4double(3*nCharged*(nCharged-1)*(nParticles-nCharged))/denominator;  
   }
   return rj;
 }
 
-////////////////////////////////////////////////////////////////////////////////////
+////////////////////////////////////////////////////////////////////////////////
 //J. M. Quesada (Dec 2007-June 2008): New inverse reaction cross sections 
 //OPT=0 Dostrovski's parameterization
@@ -87 +78 @@
 //OPT=3,4 Kalbach's parameterization 
 // 
-G4double G4PreCompoundHe3::CrossSection(const  G4double K)
-{
-  ResidualA=GetRestA();
-  ResidualZ=GetRestZ(); 
-  theA=GetA();
-  theZ=GetZ();
-  ResidualAthrd=std::pow(ResidualA,0.33333);
-  FragmentA=GetA()+GetRestA();
-  FragmentAthrd=std::pow(FragmentA,0.33333);
-
+G4double G4PreCompoundHe3::CrossSection(G4double K)
+{
+  ResidualA = GetRestA();
+  ResidualZ = GetRestZ(); 
+  theA = GetA();
+  theZ = GetZ();
+  ResidualAthrd = ResidualA13();
+  FragmentA = theA + ResidualA;
+  FragmentAthrd = g4pow->Z13(FragmentA);
 
   if (OPTxs==0) return GetOpt0( K);
@@ -109 +99 @@
 }
 
-// *********************** OPT=0 : Dostrovski's cross section  *****************************
-
-G4double G4PreCompoundHe3::GetOpt0(const  G4double K)
-{
-  const G4double r0 = G4PreCompoundParameters::GetAddress()->Getr0();
-  // cross section is now given in mb (r0 is in mm) for the sake of consistency
-  //with the rest of the options
-  return 1.e+25*pi*(r0*ResidualAthrd)*(r0*ResidualAthrd)*GetAlpha()*(1.+GetBeta()/K);
-}
-//
-//----------------
-//
 G4double G4PreCompoundHe3::GetAlpha()
 {
   G4double C = 0.0;
-  G4double aZ = GetZ() + GetRestZ();
+  G4int aZ = theZ + ResidualZ;
   if (aZ <= 30) 
     {
@@ -131 +109 @@
   else if (aZ <= 50) 
     {
-      C = 0.1 + -((aZ-50.)/20.)*0.02;
+      C = 0.1 - (aZ - 30)*0.001;
     } 
   else if (aZ < 70) 
     {
-      C = 0.08 + -((aZ-70.)/20.)*0.02;
+      C = 0.08 - (aZ - 50)*0.001;
     }
   else 
@@ -143 +121 @@
   return 1.0 + C*(4.0/3.0);
 }
-//
-//--------------------
-// 
-G4double G4PreCompoundHe3::GetBeta() 
-{
-  return -GetCoulombBarrier();
-}
-//
-//********************* OPT=1,2 : Chatterjee's cross section ************************ 
+
+//********************* OPT=1,2 : Chatterjee's cross section *****************
 //(fitting to cross section from Bechetti & Greenles OM potential)
 
 G4double G4PreCompoundHe3::GetOpt12(const  G4double K)
 {
-
-  G4double Kc=K;
+  G4double Kc = K;
 
   // JMQ xsec is set constat above limit of validity
-  if (K>50) Kc=50;
+  if (K > 50*MeV) { Kc = 50*MeV; }
 
   G4double landa ,mu ,nu ,p , Ec,q,r,ji,xs;
 
-  G4double      p0 = -3.06;
+  G4double     p0 = -3.06;
   G4double     p1 = 278.5;
   G4double     p2 = -1389.;
@@ -179 +149 @@
   p = p0 + p1/Ec + p2/(Ec*Ec);
   landa = landa0*ResidualA + landa1;
-  mu = mu0*std::pow(ResidualA,mu1);
-  nu = std::pow(ResidualA,mu1)*(nu0 + nu1*Ec + nu2*(Ec*Ec));
+
+  G4double resmu1 = g4pow->powZ(ResidualA,mu1); 
+  mu = mu0*resmu1;
+  nu = resmu1*(nu0 + nu1*Ec + nu2*(Ec*Ec));
   q = landa - nu/(Ec*Ec) - 2*p*Ec;
   r = mu + 2*nu/Ec + p*(Ec*Ec);
@@ -198 +170 @@
 //c     ** 3he from o.m. of gibson et al
 {
-
   G4double landa, mu, nu, p , signor(1.),sig;
   G4double ec,ecsq,xnulam,etest(0.),a; 
   G4double b,ecut,cut,ecut2,geom,elab;
 
-
   G4double     flow = 1.e-18;
   G4double     spill= 1.e+18;
-
 
   G4double     p0 = -2.88;
@@ -227 +196 @@
   p = p0 + p1/ec + p2/ecsq;
   landa = landa0*ResidualA + landa1;
-  a = std::pow(ResidualA,mu1);
+  a = g4pow->powZ(ResidualA,mu1);
   mu = mu0 * a;
   nu = a* (nu0+nu1*ec+nu2*ecsq);  
   xnulam = nu / landa;
-  if (xnulam > spill) xnulam=0.;
-  if (xnulam >= flow) etest = 1.2 *std::sqrt(xnulam);
+  if (xnulam > spill) { xnulam=0.; }
+  if (xnulam >= flow) { etest = 1.2 *std::sqrt(xnulam); }
   
   a = -2.*p*ec + landa - nu/ecsq;
@@ -241 +210 @@
   ecut = (ecut-a) / (p+p);
   ecut2 = ecut;
-//JMQ 290310 for avoiding unphysical increase below minimum (at ecut)
-//ecut<0 means that there is no cut with energy axis, i.e. xs is set to 0 bellow minimum
-//  if (cut < 0.) ecut2 = ecut - 2.;
-  if (cut < 0.) ecut2 = ecut;
+  //JMQ 290310 for avoiding unphysical increase below minimum (at ecut)
+  // ecut<0 means that there is no cut with energy axis, i.e. xs is set
+  // to 0 bellow minimum
+  //  if (cut < 0.) ecut2 = ecut - 2.;
+  if (cut < 0.) { ecut2 = ecut; }
   elab = K * FragmentA / ResidualA;
   sig = 0.;
   
   if (elab <= ec) { //start for E<Ec
-    if (elab > ecut2)  sig = (p*elab*elab+a*elab+b) * signor;
+    if (elab > ecut2) { sig = (p*elab*elab+a*elab+b) * signor; }
   }           //end for E<Ec
   else {           //start for E>Ec
     sig = (landa*elab+mu+nu/elab) * signor;
     geom = 0.;
-    if (xnulam < flow || elab < etest) return sig;
+    if (xnulam < flow || elab < etest) { return sig; }
     geom = std::sqrt(theA*K);
     geom = 1.23*ResidualAthrd + ra + 4.573/geom;
@@ -263 +233 @@
   
 }
-
-//   ************************** end of cross sections ******************************* 
-
-
-
-

trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/src/G4PreCompoundIon.cc

@@ -24 +24 @@
 // ********************************************************************
 //
-// $Id: G4PreCompoundIon.cc,v 1.16 2009/02/10 16:01:37 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-04-beta-01 $
+// $Id: G4PreCompoundIon.cc,v 1.17 2010/08/28 15:16:55 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03-ref-09 $
 //
 // -------------------------------------------------------------------
@@ -38 +38 @@
 // Modified:  
 // 10.02.2009 J. M. Quesada fixed bug in density level of light fragments  
+// 20.08.2010 V.Ivanchenko added G4Pow and G4PreCompoundParameters pointers
+//                         use int Z and A and cleanup
 //
 
 #include "G4PreCompoundIon.hh"
-#include "G4PreCompoundParameters.hh"
 
-G4bool G4PreCompoundIon::IsItPossible(const G4Fragment& aFragment) 
+G4PreCompoundIon::
+G4PreCompoundIon(const G4ParticleDefinition* part,
+                 G4VCoulombBarrier* aCoulombBarrier)
+  : G4PreCompoundFragment(part,aCoulombBarrier)
 {
-  G4int pplus = aFragment.GetNumberOfCharged();   
-  G4int pneut = aFragment.GetNumberOfParticles()-pplus;
-  return (pneut >= (GetA()-GetZ()) && pplus >= GetZ());
+  G4double r0 = theParameters->Getr0();
+  fact = 0.75*CLHEP::millibarn/(CLHEP::pi*r0*r0*r0);
 }
 
+G4PreCompoundIon::~G4PreCompoundIon()
+{}
+
 G4double G4PreCompoundIon::
-ProbabilityDistributionFunction(const G4double eKin, 
+ProbabilityDistributionFunction(G4double eKin, 
                                 const G4Fragment& aFragment)
 {
-  if ( !IsItPossible(aFragment) ) return 0.0;
-  
-  const G4double r0 = G4PreCompoundParameters::GetAddress()->Getr0();
+  if ( !IsItPossible(aFragment) ) { return 0.0; }
+  G4double efinal = eKin + GetBindingEnergy();
+  //G4cout << "Efinal= " << efinal << " Ekin= " << eKin << G4endl;
+  if(efinal <= 0.0 ) { return 0.0; } 
 
   G4double U = aFragment.GetExcitationEnergy();
-  G4double P = aFragment.GetNumberOfParticles();
-  G4double H = aFragment.GetNumberOfHoles();
-  G4double N = P + H;
+  G4int P = aFragment.GetNumberOfParticles();
+  G4int H = aFragment.GetNumberOfHoles();
+  G4int A = GetA();
+  G4int N = P + H;
 
-  G4double g0 = (6.0/pi2)*aFragment.GetA() * 
-    G4PreCompoundParameters::GetAddress()->GetLevelDensity();
- 
-  G4double g1 = (6.0/pi2)*GetRestA() * 
-    G4PreCompoundParameters::GetAddress()->GetLevelDensity();
+  G4double g0 = (6.0/pi2)*aFragment.GetA_asInt()*theParameters->GetLevelDensity();
+  G4double g1 = (6.0/pi2)*GetRestA()*theParameters->GetLevelDensity();
 
   //JMQ 06/02/209  This is  THE BUG that was killing cluster emission
@@ -75 +80 @@
   G4double gj = g1;
 
-  G4double A0 = ((P*P+H*H+P-H)/4.0 - H/2.0)/g0; 
+  G4double A0 = G4double(P*P+H*H+P-3*H)/(4.0*g0);
+  G4double A1 = std::max(0.0,(A0*g0 + A*(A-2*P-1)*0.25)/g1); 
 
-  G4double A1 = std::max(0.0,(A0*g0 + GetA()*(GetA()-2.0*P-1.0)/4.0)/g1); 
-
-  G4double Aj = GetA()*(GetA()+1.0)/4.0/gj; 
-
-
-  G4double E0 = std::max(0.0,U - A0);
-  if (E0 == 0.0) return 0.0;
+  G4double E0 = U - A0;
+  //G4cout << "E0= " << E0 << G4endl;
+  if (E0 <= 0.0) { return 0.0; }
 
   G4double E1 = (std::max(0.0,GetMaximalKineticEnergy() - eKin - A1)); 
 
-  G4double Ej = std::max(0.0,eKin + GetBindingEnergy() -Aj); 
+  G4double Aj = A*(A+1)/(4.0*gj); 
+  G4double Ej = std::max(0.0,efinal - Aj); 
+
+  G4double rj = GetRj(P, aFragment.GetNumberOfCharged());
+  G4double xs = CrossSection(eKin);
+
+  //G4cout << "rj= " << rj << " xs= " << xs << G4endl;
 
   // JMQ 10/02/09 reshaping of the formula (unnecessary std::pow elimitated)
+  /*
+  G4double r0 = theParameters->Getr0();
   G4double pA = (3.0/4.0) * std::sqrt(std::max(0.0, 2.0/(GetReducedMass()*
-                (eKin+GetBindingEnergy()))))/(pi * r0 * r0 *r0* GetRestA())* 
+                (eKin+GetBindingEnergy()))))/(pi * r0 * r0 *r0* GetRestA())* 
                 eKin*CrossSection(eKin)*millibarn* 
-                CoalescenceFactor(aFragment.GetA()) * FactorialFactor(N,P)* 
-                GetRj(aFragment.GetNumberOfParticles(), aFragment.GetNumberOfCharged());
+                CoalescenceFactor(aFragment.GetA_asInt()) * FactorialFactor(N,P)* 
+       GetRj(aFragment.GetNumberOfParticles(), aFragment.GetNumberOfCharged());
 
   G4double pB = std::pow((g1*E1)/(g0*E0),N-GetA()-1.0)*(g1/g0);
- 
   G4double pC = std::pow((gj*Ej)/(g0*E0),GetA()-1.0)*(gj/g0)/E0; 
-
-  G4double Probability = pA * pB * pC;
-
-  return Probability;
+  pA *= pB * pC;
+  */
+  
+  G4double pA = fact*eKin*xs*rj 
+    * CoalescenceFactor(aFragment.GetA_asInt()) * FactorialFactor(N,P)
+    * std::sqrt(2.0/(GetReducedMass()*efinal)) 
+    * g4pow->powN(g1*E1/(g0*E0), N-A-1)
+    * g4pow->powN(gj*Ej/(g0*E0), A-1)*gj*g1/(g0*g0*E0*GetRestA()); 
+   
+  return pA;
 }

trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/src/G4PreCompoundModel.cc

@@ -24 +24 @@
 // ********************************************************************
 //
-//
-// $Id: G4PreCompoundModel.cc,v 1.20 2010/06/11 17:26:43 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-04-beta-01 $
+// $Id: G4PreCompoundModel.cc,v 1.25 2010/10/11 13:54:59 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03-ref-09 $
 //
 // by V. Lara
 //
-//J. M. Quesada (Apr.08). Several changes. Soft cut-off switched off. 
-//(May. 08). Protection against non-physical preeq. transitional regime has 
-// been set
-//
-// Modif (03 September 2008) by J. M. Quesada for external choice of inverse 
-// cross section option
-// JMQ (06 September 2008) Also external choices have been added for:
+// Modified:
+// 01.04.2008 J.M.Quesada Several changes. Soft cut-off switched off. 
+// 01.05.2008 J.M.Quesada Protection against non-physical preeq. 
+//                        transitional regime has been set
+// 03.09.2008 J.M.Quesada for external choice of inverse cross section option
+// 06.09.2008 J.M.Quesada Also external choices have been added for:
 //                      - superimposed Coulomb barrier (useSICB=true) 
 //                      - "never go back"  hipothesis (useNGB=true) 
 //                      - soft cutoff from preeq. to equlibrium (useSCO=true)
 //                      - CEM transition probabilities (useCEMtr=true)  
-
+// 20.08.2010 V.Ivanchenko Cleanup of the code: 
+//                      - integer Z and A;
+//                      - emission and transition classes created at initialisation
+//                      - options are set at initialisation
+//                      - do not use copy-constructors for G4Fragment  
 
 #include "G4PreCompoundModel.hh"
@@ -48 +50 @@
 #include "G4GNASHTransitions.hh"
 #include "G4ParticleDefinition.hh"
-
+#include "G4Proton.hh"
+#include "G4Neutron.hh"
+
+#include "G4NucleiProperties.hh"
+#include "G4PreCompoundParameters.hh"
+#include "Randomize.hh"
+#include "G4DynamicParticle.hh"
+#include "G4ParticleTypes.hh"
+#include "G4ParticleTable.hh"
+#include "G4LorentzVector.hh"
 
 #ifdef PRECOMPOUND_TEST
@@ -58 +69 @@
   : G4VPreCompoundModel(value), useHETCEmission(false), useGNASHTransition(false), 
     OPTxs(3), useSICB(false), useNGB(false), useSCO(false), useCEMtr(true) 
-{}
+{
+  theParameters = G4PreCompoundParameters::GetAddress();
+
+  theEmission = new G4PreCompoundEmission();
+  if(useHETCEmission) { theEmission->SetHETCModel(); }
+  else { theEmission->SetDefaultModel(); }
+  theEmission->SetOPTxs(OPTxs);
+  theEmission->UseSICB(useSICB);
+
+  if(useGNASHTransition) { theTransition = new G4GNASHTransitions; }
+  else { theTransition = new G4PreCompoundTransitions(); }
+  theTransition->UseNGB(useNGB);
+  theTransition->UseCEMtr(useCEMtr);
+
+  proton = G4Proton::Proton();
+  neutron = G4Neutron::Neutron();
+}
 
 G4PreCompoundModel::~G4PreCompoundModel() 
-{}
-
-G4PreCompoundModel::G4PreCompoundModel() 
-  : G4VPreCompoundModel(0), useHETCEmission(false), useGNASHTransition(false), 
-    OPTxs(3), useSICB(false), useNGB(false), useSCO(false), useCEMtr(true) 
-{}
-
-G4PreCompoundModel::G4PreCompoundModel(const G4PreCompoundModel &) 
-: G4VPreCompoundModel() 
-{}
-
-const G4PreCompoundModel & G4PreCompoundModel::operator=(const G4PreCompoundModel &)
 {
-  throw G4HadronicException(__FILE__, __LINE__, "G4PreCompoundModel::operator= meant to not be accessable");
-  return *this;
-}
-
-
-G4bool G4PreCompoundModel::operator==(const G4PreCompoundModel &) const
-{
-  return false;
-}
-
-G4bool G4PreCompoundModel::operator!=(const G4PreCompoundModel &) const
-{
-  return true;
-}
-
-
-
-// Additional Declarations
-
-G4HadFinalState * G4PreCompoundModel::ApplyYourself(const G4HadProjectile & thePrimary,
-                                                    G4Nucleus & theNucleus)
+  delete theEmission;
+  delete theTransition;
+}
+
+/////////////////////////////////////////////////////////////////////////////////////////
+
+G4HadFinalState* G4PreCompoundModel::ApplyYourself(const G4HadProjectile & thePrimary,
+                                                   G4Nucleus & theNucleus)
 {  
+  const G4ParticleDefinition* primary = thePrimary.GetDefinition();
+  if(primary != neutron && primary != proton) {
+    std::ostringstream errOs;
+    errOs << "BAD primary type in G4PreCompoundModel: " 
+          << primary->GetParticleName() <<G4endl;
+    throw G4HadronicException(__FILE__, __LINE__, errOs.str());
+  }
+  G4int Zp = 0;
+  G4int Ap = 1;
+  if(primary == proton) { Zp = 1; }
+
+  G4int A = theNucleus.GetA_asInt();
+  G4int Z = theNucleus.GetZ_asInt();
+
+  // 4-Momentum
+  G4LorentzVector p = thePrimary.Get4Momentum();
+  G4double mass = G4NucleiProperties::GetNuclearMass(A, Z);
+  p += G4LorentzVector(0.0,0.0,0.0,mass);
+
   // prepare fragment
-  G4Fragment anInitialState;
-  // This si for GNASH transitions
-  anInitialState.SetParticleDefinition(const_cast<G4ParticleDefinition *>(thePrimary.GetDefinition()));
-
-  G4int anA=static_cast<G4int>(theNucleus.GetN());
-  anA += thePrimary.GetDefinition()->GetBaryonNumber();
-
-  anInitialState.SetA(anA);
-  
-  G4int aZ=static_cast<G4int>(theNucleus.GetZ());
-  aZ += static_cast<G4int>(thePrimary.GetDefinition()->GetPDGCharge());
-
-  anInitialState.SetZ(aZ);
-  
-  // Assume the projectile is a nucleon
-  
+  G4Fragment anInitialState(A + Ap, Z + Zp, p);
+  anInitialState.SetNumberOfParticles(2);
+  anInitialState.SetNumberOfHoles(1);
+  anInitialState.SetNumberOfCharged(1);
+  anInitialState.SetCreationTime(thePrimary.GetGlobalTime());
+  /*
   // Number of Excited Particles
   anInitialState.SetNumberOfParticles(1+thePrimary.GetDefinition()->GetBaryonNumber());
@@ -129 +140 @@
   // Number of Holes 
   anInitialState.SetNumberOfHoles(1);
-  
-  // pre-compound nucleus energy.
-  G4double anEnergy = 0;
-  G4double nucleusMass =  G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(static_cast<G4int>(theNucleus.GetZ()),
-                                                                                         static_cast<G4int>(theNucleus.GetN()));
-  anEnergy =  nucleusMass + thePrimary.GetTotalEnergy();
-  
-  // Momentum
-  G4ThreeVector p = thePrimary.Get4Momentum().vect();
-
-  // 4-momentum
-  G4LorentzVector momentum(p, anEnergy);
-  anInitialState.SetMomentum(momentum);
-  
-#ifdef PRECOMPOUND_TEST
-  G4PreCompoundModel::theInitialFragmentForTest = anInitialState;
-#endif
+  */
   
   // call excitation handler
-  const G4Fragment aFragment(anInitialState);
-  G4ReactionProductVector * result = DeExcite(aFragment);
+  //  const G4Fragment aFragment(anInitialState);
+  G4ReactionProductVector * result = DeExcite(anInitialState);
 
 #ifdef PRECOMPOUND_TEST
@@ -180 +175 @@
 }
 
-
 /////////////////////////////////////////////////////////////////////////////////////////
-/////////////////////////////////////////////////////////////////////////////////////////
-
-G4ReactionProductVector* G4PreCompoundModel::DeExcite(const G4Fragment & theInitialState) const
+
+G4ReactionProductVector* G4PreCompoundModel::DeExcite(G4Fragment& aFragment)
 {
-  
   G4ReactionProductVector * Result = new G4ReactionProductVector;
-  
-  // Copy of the initial state 
-  G4Fragment aFragment(theInitialState);
-
-  if (aFragment.GetExcitationEnergy() < 10*eV)
-    {
-      // Perform Equilibrium Emission
-      PerformEquilibriumEmission(aFragment,Result);
-      return Result;
-    }
-  
-  if (aFragment.GetA() < 5) {
-    G4ReactionProduct * theRP = new G4ReactionProduct(G4ParticleTable::GetParticleTable()->
-                                                      GetIon(static_cast<G4int>(aFragment.GetZ()),
-                                                             static_cast<G4int>(aFragment.GetA()),
-                                                             aFragment.GetExcitationEnergy()));
-    theRP->SetMomentum(aFragment.GetMomentum().vect());
-    theRP->SetTotalEnergy(aFragment.GetMomentum().e());   
-    Result->push_back(theRP);
+  G4double Eex = aFragment.GetExcitationEnergy();
+  G4int A = aFragment.GetA_asInt(); 
+
+  //G4cout << "### G4PreCompoundModel::DeExcite" << G4endl;
+  //G4cout << aFragment << G4endl;
+ 
+  // Perform Equilibrium Emission 
+  if (A < 5 || Eex < keV /*|| Eex > 3.*MeV*A*/) {
+    PerformEquilibriumEmission(aFragment, Result);
     return Result;
   }
   
-  G4PreCompoundEmission aEmission;
-  if (useHETCEmission) aEmission.SetHETCModel();
-  aEmission.SetUp(theInitialState);
-  
-  //for cross section options
-  
-  if (OPTxs!= 0 && OPTxs!=1 && OPTxs !=2 && OPTxs !=3 && OPTxs !=4  ) {
-    std::ostringstream errOs;
-    errOs << "BAD CROSS SECTION OPTION in G4PreCompoundModel.cc !!"  <<G4endl;
-    throw G4HadronicException(__FILE__, __LINE__, errOs.str());}
-  else aEmission.SetOPTxs(OPTxs); 
-  
-  //for the choice of superimposed Coulomb Barrier for inverse cross sections
-  
-   aEmission.UseSICB(useSICB);
- 
-  
-  //----------
-  
-  G4VPreCompoundTransitions * aTransition = 0;
-  if (useGNASHTransition) 
-    {
-      aTransition = new G4GNASHTransitions;
-    }
-  else 
-    {
-      aTransition = new G4PreCompoundTransitions;
-      // for the choice of "never go back" hypothesis and CEM transition probabilities  
-      if (useNGB) aTransition->UseNGB(useNGB);
-      if (useCEMtr) aTransition->UseCEMtr(useCEMtr);
-    }
-  
-  // Main loop. It is performed until equilibrium deexcitation.
-  //G4int fragment=0;
-  
+  // main loop  
   for (;;) {
     
@@ -252 +200 @@
     
     // Initialize fragment according with the nucleus parameters
-    aEmission.Initialize(aFragment);
-    
-    
-    
-    G4double g = (6.0/pi2)*aFragment.GetA()*
-      G4PreCompoundParameters::GetAddress()->GetLevelDensity();
-    
-    
-    
-    
-    G4int EquilibriumExcitonNumber = static_cast<G4int>(std::sqrt(2.0*g*aFragment.GetExcitationEnergy())+ 0.5);
-//   
-//    G4cout<<"Neq="<<EquilibriumExcitonNumber<<G4endl;
-//
-// J. M. Quesada (Jan. 08)  equilibrium hole number could be used as preeq.- evap. delimiter (IAEA report)
-//    G4int EquilibriumHoleNumber = static_cast<G4int>(0.2*std::sqrt(g*aFragment.GetExcitationEnergy())+ 0.5);
-    
-// Loop for transitions, it is performed while there are preequilibrium transitions.
+    theEmission->Initialize(aFragment);
+    
+    G4double g = (6.0/pi2)*aFragment.GetA_asInt()*theParameters->GetLevelDensity();
+    
+    G4int EquilibriumExcitonNumber = 
+      static_cast<G4int>(std::sqrt(2.0*g*aFragment.GetExcitationEnergy())+ 0.5);
+    //   
+    //    G4cout<<"Neq="<<EquilibriumExcitonNumber<<G4endl;
+    //
+    // J. M. Quesada (Jan. 08)  equilibrium hole number could be used as preeq.
+    // evap. delimiter (IAEA report)
+    
+    // Loop for transitions, it is performed while there are preequilibrium transitions.
     G4bool ThereIsTransition = false;
     
@@ -278 +221 @@
     //        G4cout<<" Ex. Energy="<<aFragment.GetExcitationEnergy()<<G4endl;
     //        G4cout<<"N. excitons="<<NE<<"  N. Part="<<NP<<"N. Holes ="<<NH<<G4endl;
-        
-    
     //G4int transition=0;
-    do 
-      {
-        //transition++;
-        //G4cout<<"transition number .."<<transition<<G4endl;
-        //G4cout<<" n ="<<aFragment.GetNumberOfExcitons()<<G4endl;
-        G4bool go_ahead = false;
-        // soft cutoff criterium as an "ad-hoc" solution to force increase in  evaporation  
-        //       G4double test = static_cast<G4double>(aFragment.GetNumberOfHoles());
-        G4double test = static_cast<G4double>(aFragment.GetNumberOfExcitons());
-        
-        
-        if (test < EquilibriumExcitonNumber) go_ahead=true;
-        //J. M. Quesada (Apr. 08): soft-cutoff switched off by default
-        if (useSCO) {
-          if (test < EquilibriumExcitonNumber)
-            //  if (test < EquilibriumHoleNumber)
-            {
-              test /= static_cast<G4double>(EquilibriumExcitonNumber); 
-              //     test /= static_cast<G4double>(EquilibriumHoleNumber);
-              test -= 1.0;
+    do {
+      //transition++;
+      //G4cout<<"transition number .."<<transition<<G4endl;
+      //G4cout<<" n ="<<aFragment.GetNumberOfExcitons()<<G4endl;
+      G4bool go_ahead = false;
+      // soft cutoff criterium as an "ad-hoc" solution to force increase in  evaporation  
+      //       G4double test = static_cast<G4double>(aFragment.GetNumberOfHoles());
+      G4int test = aFragment.GetNumberOfExcitons();
+      if (test < EquilibriumExcitonNumber) { go_ahead=true; }
+
+      //J. M. Quesada (Apr. 08): soft-cutoff switched off by default
+      if (useSCO) {
+        if (test < EquilibriumExcitonNumber)
+          {
+            G4double x = G4double(test)/G4double(EquilibriumExcitonNumber) - 1; 
+            if( G4UniformRand() < 1.0 -  std::exp(-x*x/0.32) ) { go_ahead = true; }
+            /*
               test = test*test;
               test /= 0.32;
               test = 1.0 - std::exp(-test);
               go_ahead = (G4UniformRand() < test);
-              
-            }
-        } 
+            */
+          }
+      } 
         
-        //JMQ: WARNING:  CalculateProbability MUST be called prior to Get methods !! (O values would be returned otherwise)
-        G4double TotalTransitionProbability = aTransition->CalculateProbability(aFragment);
-        G4double P1=aTransition->GetTransitionProb1();
-        G4double P2=aTransition->GetTransitionProb2();
-        G4double P3=aTransition->GetTransitionProb3();
-        //       G4cout<<"P1="<<P1<<" P2="<<P2<<"  P3="<<P3<<G4endl;
+      // JMQ: WARNING:  CalculateProbability MUST be called prior to Get methods !! 
+      // (O values would be returned otherwise)
+      G4double TotalTransitionProbability = 
+        theTransition->CalculateProbability(aFragment);
+      G4double P1 = theTransition->GetTransitionProb1();
+      G4double P2 = theTransition->GetTransitionProb2();
+      G4double P3 = theTransition->GetTransitionProb3();
+      //G4cout<<"P1="<<P1<<" P2="<<P2<<"  P3="<<P3<<G4endl;
+      
+      //J.M. Quesada (May. 08). Physical criterium (lamdas)  PREVAILS over 
+      //                        approximation (critical exciton number)
+      if(P1 <= P2+P3) { go_ahead = false; }
         
-        
-        //J.M. Quesada (May. 08). Physical criterium (lamdas)  PREVAILS over approximation (critical exciton number)
-        if(P1<=(P2+P3)) go_ahead=false;
-        
-        if (go_ahead &&  aFragment.GetA() > 4)
-          {                             
-            G4double TotalEmissionProbability = aEmission.GetTotalProbability(aFragment);
-            //
-            //  G4cout<<"TotalEmissionProbability="<<TotalEmissionProbability<<G4endl;
-            //
-            // Check if number of excitons is greater than 0
-            // else perform equilibrium emission
-            if (aFragment.GetNumberOfExcitons() <= 0) 
-              {
-                // Perform Equilibrium Emission
-#ifdef debug // ------------- debug -----------------------------------------
-                CheckConservation(theInitialState,aFragment,Result);
-#endif // ------------------- debug -----------------------------------------
-                PerformEquilibriumEmission(aFragment,Result);
-                delete aTransition;
-                return Result;
-              }
+      if (go_ahead &&  aFragment.GetA_asInt() > 4) 
+        {
+                                
+          G4double TotalEmissionProbability = 
+            theEmission->GetTotalProbability(aFragment);
+          //
+          //  G4cout<<"TotalEmissionProbability="<<TotalEmissionProbability<<G4endl;
+          //
+          // Check if number of excitons is greater than 0
+          // else perform equilibrium emission
+          if (aFragment.GetNumberOfExcitons() <= 0) 
+            {
+              PerformEquilibriumEmission(aFragment,Result);
+              return Result;
+            }
             
-            //      G4PreCompoundTransitions aTransition(aFragment);
+          //J.M.Quesada (May 08) this has already been done in order to decide  
+          //                     what to do (preeq-eq) 
+          // Sum of all probabilities
+          G4double TotalProbability = TotalEmissionProbability 
+            + TotalTransitionProbability;
             
-            //J.M.Quesada (May 08) this has already been done in order to decide what to do (preeq-eq) 
-            // Sum of transition probabilities
-            //      G4double TotalTransitionProbability = aTransition->CalculateProbability(aFragment);
-            
-            // Sum of all probabilities
-            G4double TotalProbability = TotalEmissionProbability + TotalTransitionProbability;
-            
-            // Select subprocess
-            if (G4UniformRand() > TotalEmissionProbability/TotalProbability) 
-              {
-                // It will be transition to state with a new number of excitons
-                ThereIsTransition = true;               
-                // Perform the transition
-                aFragment = aTransition->PerformTransition(aFragment);
-              } 
-            else 
-              {
-                // It will be fragment emission
-                ThereIsTransition = false;
-                Result->push_back(aEmission.PerformEmission(aFragment));
-              }
-          } 
-        else 
-          {
-            // Perform Equilibrium Emission
-#ifdef debug
-            CheckConservation(theInitialState,aFragment,Result);
-#endif
-            PerformEquilibriumEmission(aFragment,Result);
-            delete aTransition;
-            return Result;
-          }
-      } while (ThereIsTransition);   // end of do loop
+          // Select subprocess
+          if (TotalProbability*G4UniformRand() > TotalEmissionProbability) 
+            {
+              // It will be transition to state with a new number of excitons
+              ThereIsTransition = true;         
+              // Perform the transition
+              theTransition->PerformTransition(aFragment);
+            } 
+          else 
+            {
+              // It will be fragment emission
+              ThereIsTransition = false;
+              Result->push_back(theEmission->PerformEmission(aFragment));
+            }
+        } 
+      else 
+        {
+          PerformEquilibriumEmission(aFragment,Result);
+          return Result;
+        }
+    } while (ThereIsTransition);   // end of do loop
   } // end of for (;;) loop
-}
-
-
-
-
-void G4PreCompoundModel::PerformEquilibriumEmission(const G4Fragment & aFragment,
-                                                    G4ReactionProductVector * Result) const 
-{
-  G4ReactionProductVector * theEquilibriumResult;
-
-  theEquilibriumResult = GetExcitationHandler()->BreakItUp(aFragment);
-  
-  Result->insert(Result->end(),theEquilibriumResult->begin(), theEquilibriumResult->end());
-
-  delete theEquilibriumResult;
-  return;
-}
-
+  return Result;
+}
+
+/////////////////////////////////////////////////////////////////////////////////////////
+//       Initialisation
+/////////////////////////////////////////////////////////////////////////////////////////
+
+void G4PreCompoundModel::UseHETCEmission() 
+{ 
+  useHETCEmission = true; 
+  theEmission->SetHETCModel();
+}
+
+void G4PreCompoundModel::UseDefaultEmission() 
+{ 
+  useHETCEmission = false; 
+  theEmission->SetDefaultModel();
+}
+
+void G4PreCompoundModel::UseGNASHTransition() { 
+  useGNASHTransition = true; 
+  delete theTransition;
+  theTransition = new G4GNASHTransitions;
+  theTransition->UseNGB(useNGB);
+  theTransition->UseCEMtr(useCEMtr);
+}
+
+void G4PreCompoundModel::UseDefaultTransition() { 
+  useGNASHTransition = false; 
+  delete theTransition;
+  theTransition = new G4PreCompoundTransitions();
+  theTransition->UseNGB(useNGB);
+  theTransition->UseCEMtr(useCEMtr);
+}
+
+void G4PreCompoundModel::SetOPTxs(G4int opt) 
+{ 
+  OPTxs = opt; 
+  theEmission->SetOPTxs(OPTxs);
+}
+
+void G4PreCompoundModel::UseSICB() 
+{ 
+  useSICB = true; 
+  theEmission->UseSICB(useSICB);
+}
+
+void G4PreCompoundModel::UseNGB()  
+{ 
+  useNGB = true; 
+}
+
+void G4PreCompoundModel::UseSCO()  
+{ 
+  useSCO = true; 
+}
+
+void G4PreCompoundModel::UseCEMtr() 
+{ 
+  useCEMtr = true; 
+}
+
+/////////////////////////////////////////////////////////////////////////////////////////
 
 #ifdef debug

trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/src/G4PreCompoundNeutron.cc

@@ -24 +24 @@
 // ********************************************************************
 //
-//
-// $Id: G4PreCompoundNeutron.cc,v 1.4 2009/02/11 18:06:00 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-04-beta-01 $
+// $Id: G4PreCompoundNeutron.cc,v 1.5 2010/08/28 15:16:55 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03-ref-09 $
 //
 // -------------------------------------------------------------------
@@ -40 +39 @@
 // 21.08.2008 J. M. Quesada add choice of options  
 // 10.02.2009 J. M. Quesada set default opt3
+// 20.08.2010 V.Ivanchenko added G4Pow and G4PreCompoundParameters pointers
+//                         use int Z and A and cleanup
 // 
 
 #include "G4PreCompoundNeutron.hh"
-
-G4ReactionProduct * G4PreCompoundNeutron::GetReactionProduct() const
-{
-  G4ReactionProduct * theReactionProduct = 
-    new G4ReactionProduct(G4Neutron::NeutronDefinition());
-  theReactionProduct->SetMomentum(GetMomentum().vect());
-  theReactionProduct->SetTotalEnergy(GetMomentum().e());
-#ifdef PRECOMPOUND_TEST
-  theReactionProduct->SetCreatorModel("G4PrecompoundModel");
-#endif
-  return theReactionProduct;
-}
-
-G4double G4PreCompoundNeutron::GetRj(const G4int NumberParticles, const G4int NumberCharged)
+#include "G4Neutron.hh"
+
+G4PreCompoundNeutron::G4PreCompoundNeutron()
+  : G4PreCompoundNucleon(G4Neutron::Neutron(), &theNeutronCoulombBarrier)
+{}
+
+G4PreCompoundNeutron::~G4PreCompoundNeutron()
+{}
+
+G4double G4PreCompoundNeutron::GetRj(G4int nParticles, G4int nCharged)
 {
   G4double rj = 0.0;
-  if(NumberParticles > 0) rj = static_cast<G4double>(NumberParticles - NumberCharged)/
-                            static_cast<G4double>(NumberParticles);
+  if(nParticles > 0) { 
+    rj = static_cast<G4double>(nParticles - nCharged)/
+      static_cast<G4double>(nParticles);
+  }
   return rj;
 }
@@ -72 +71 @@
 G4double G4PreCompoundNeutron::CrossSection(const  G4double K)
 {
-  ResidualA=GetRestA();
-  ResidualZ=GetRestZ(); 
-  theA=GetA();
-  theZ=GetZ();
-  ResidualAthrd=std::pow(ResidualA,0.33333);
-  FragmentA=GetA()+GetRestA();
-  FragmentAthrd=std::pow(FragmentA,0.33333);
-
-  if (OPTxs==0) return GetOpt0( K);
-  else if( OPTxs==1 || OPTxs==2) return GetOpt12( K);
-  else if (OPTxs==3 || OPTxs==4)  return GetOpt34( K);
+  ResidualA = GetRestA();
+  ResidualZ = GetRestZ(); 
+  theA = GetA();
+  theZ = GetZ();
+  ResidualAthrd = ResidualA13();
+  FragmentA = theA + ResidualA;
+  FragmentAthrd = g4pow->Z13(FragmentA);
+
+  if (OPTxs==0) { return GetOpt0( K); }
+  else if( OPTxs==1 || OPTxs==2) { return GetOpt12( K); }
+  else if (OPTxs==3 || OPTxs==4) { return GetOpt34( K); }
   else{
     std::ostringstream errOs;
@@ -91 +90 @@
 }
 
-// *********************** OPT=0 : Dostrovski's cross section  *****************************
-
-G4double G4PreCompoundNeutron::GetOpt0(const  G4double K)
-{
-  
-  const G4double r0 = G4PreCompoundParameters::GetAddress()->Getr0();
-  // cross section is now given in mb (r0 is in mm) for the sake of consistency
-  //with the rest of the options
-  return 1.e+25*pi*(r0*ResidualAthrd)*(r0*ResidualAthrd)*GetAlpha()*(1.+GetBeta()/K);
-}
-//
-//-------
-//
 G4double G4PreCompoundNeutron::GetAlpha()
 {
-  //   return 0.76+2.2/std::pow(GetRestA(),1.0/3.0);
   return 0.76+2.2/ResidualAthrd;
 }
-//
-//------------
-//
+
 G4double G4PreCompoundNeutron::GetBeta() 
 {
@@ -117 +100 @@
   return (2.12/(ResidualAthrd*ResidualAthrd)-0.05)*MeV/GetAlpha();
 }
-//
-
-//********************* OPT=1,2 : Chatterjee's cross section ************************ 
+
+//********************* OPT=1,2 : Chatterjee's cross section *******************
 //(fitting to cross section from Bechetti & Greenles OM potential)
 
-G4double G4PreCompoundNeutron::GetOpt12(const  G4double K)
-{
-
+G4double G4PreCompoundNeutron::GetOpt12(G4double K)
+{
   G4double Kc=K;
 
-// Pramana (Bechetti & Greenles) for neutrons is chosen 
-
-// JMQ  xsec is set constat above limit of validity
- if (K>50) Kc=50;
+  // Pramana (Bechetti & Greenles) for neutrons is chosen 
+
+  // JMQ  xsec is set constat above limit of validity
+  if (K > 50*MeV) { Kc = 50*MeV; }
 
  G4double landa, landa0, landa1, mu, mu0, mu1,nu, nu0, nu1, nu2,xs;
@@ -155 +136 @@
 }
 
-
 // *********** OPT=3,4 : Kalbach's cross sections (from PRECO code)*************
-G4double G4PreCompoundNeutron::GetOpt34(const  G4double K)
-{
-  
+G4double G4PreCompoundNeutron::GetOpt34(G4double K)
+{
   G4double landa, landa0, landa1, mu, mu0, mu1,nu, nu0, nu1, nu2;
   G4double p, p0, p1, p2;
@@ -177 +156 @@
   p2=0.;
 
-
   flow = 1.e-18;
   spill= 1.0e+18; 
 
   // PRECO xs for neutrons is choosen
-
   p0 = -312.;
   p1= 0.;
@@ -194 +171 @@
   nu2 = 1280.8; 
 
-  if (ResidualA < 40.) signor=0.7+ResidualA*0.0075;
-  if (ResidualA > 210.) signor = 1. + (ResidualA-210.)/250.;
+  if (ResidualA < 40)  { signor =0.7 + ResidualA*0.0075; }
+  if (ResidualA > 210) { signor = 1. + (ResidualA-210)/250.; }
   landa = landa0/ResidualAthrd + landa1;
   mu = mu0*ResidualAthrd + mu1*ResidualAthrd*ResidualAthrd;
@@ -201 +178 @@
 
   // JMQ very low energy behaviour corrected (problem  for A (apprx.)>60)
-  if (nu < 0.)nu=-nu;
+  if (nu < 0.) { nu=-nu; }
 
   ec = 0.5;
@@ -216 +193 @@
   ecut = 0.;
   cut = a*a - 4.*p*b;
-  if (cut > 0.) ecut = std::sqrt(cut);
+  if (cut > 0.) { ecut = std::sqrt(cut); }
   ecut = (ecut-a) / (p+p);
   ecut2 = ecut;
-  if (cut < 0.) ecut2 = ecut - 2.;
-  elab = K * FragmentA / ResidualA;
+  if (cut < 0.) { ecut2 = ecut - 2.; }
+  elab = K * FragmentA / G4double(ResidualA);
   sig = 0.;
   if (elab <= ec) { //start for E<Ec 
-    if (elab > ecut2)  sig = (p*elab*elab+a*elab+b) * signor;                
+    if (elab > ecut2) { sig = (p*elab*elab+a*elab+b) * signor; }                
   }              //end for E<Ec
   else {           //start for  E>Ec
     sig = (landa*elab+mu+nu/elab) * signor;
     geom = 0.;
-    if (xnulam < flow || elab < etest) return sig;
+    if (xnulam < flow || elab < etest) { return sig; }
     geom = std::sqrt(theA*K);
     geom = 1.23*ResidualAthrd + ra + 4.573/geom;
@@ -235 +212 @@
 
   }
-  return sig;}
-
-//   ************************** end of cross sections ******************************* 
-
-
+  return sig;
+}

trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/src/G4PreCompoundNucleon.cc

@@ -24 +24 @@
 // ********************************************************************
 //
-//
-// $Id: G4PreCompoundNucleon.cc,v 1.13 2009/02/11 18:06:00 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-04-beta-01 $
+// $Id: G4PreCompoundNucleon.cc,v 1.14 2010/08/28 15:16:55 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03-ref-09 $
 //
 // -------------------------------------------------------------------
@@ -39 +38 @@
 // Modified:  
 // 10.02.2009 J. M. Quesada cleanup
+// 20.08.2010 V.Ivanchenko added G4Pow and G4PreCompoundParameters pointers
+//                         use int Z and A and cleanup
 //
 
 #include "G4PreCompoundNucleon.hh"
-#include "G4PreCompoundParameters.hh"
 
-G4bool G4PreCompoundNucleon::IsItPossible(const G4Fragment& aFragment) 
+G4PreCompoundNucleon::
+G4PreCompoundNucleon(const G4ParticleDefinition* part,
+                      G4VCoulombBarrier* aCoulombBarrier)
+  : G4PreCompoundFragment(part,aCoulombBarrier)
 {
-  G4int pplus = aFragment.GetNumberOfCharged();   
-  G4int pneut = aFragment.GetNumberOfParticles()-pplus;
-  return (pneut >= (GetA()-GetZ()) && pplus >= GetZ());
+  fact = 2*CLHEP::millibarn/(CLHEP::pi2*CLHEP::hbarc*CLHEP::hbarc*CLHEP::hbarc);
 }
 
+G4PreCompoundNucleon::~G4PreCompoundNucleon()
+{}
+
 G4double G4PreCompoundNucleon::
-ProbabilityDistributionFunction(const G4double eKin, 
+ProbabilityDistributionFunction(G4double eKin, 
                                 const G4Fragment& aFragment)
 {
-  if ( !IsItPossible(aFragment) ) return 0.0;
+  if ( !IsItPossible(aFragment) ) { return 0.0; }
 
   G4double U = aFragment.GetExcitationEnergy();
-  G4double P = aFragment.GetNumberOfParticles();
-  G4double H = aFragment.GetNumberOfHoles();
-  G4double N = P + H;
+  G4int P = aFragment.GetNumberOfParticles();
+  G4int H = aFragment.GetNumberOfHoles();
+  G4int N = P + H;
+
+  G4double g0 = (6.0/pi2)*aFragment.GetA_asInt()*theParameters->GetLevelDensity();
+  G4double g1 = (6.0/pi2)*GetRestA()*theParameters->GetLevelDensity();
   
-  G4double g0 = (6.0/pi2)*aFragment.GetA() * 
-    G4PreCompoundParameters::GetAddress()->GetLevelDensity();
- 
-  G4double g1 = (6.0/pi2)*GetRestA() * 
-    G4PreCompoundParameters::GetAddress()->GetLevelDensity();
+  G4double A0 = G4double(P*P+H*H+P-3*H)/(4.0*g0);
+  G4double A1 = (A0 - 0.5*P)/g1;  
 
-  G4double A0 = ((P*P+H*H+P-H)/4.0 - H/2.0)/g0;
+  G4double E0 = U - A0;
+  if (E0 <= 0.0) { return 0.0; }
 
-  G4double A1 = (A0 - P/2.0)/g1;
+  G4double E1 = U - eKin - GetBindingEnergy() - A1;
+  if (E1 <= 0.0) { return 0.0; }
+
+  G4double rj = GetRj(P, aFragment.GetNumberOfCharged());
+  G4double xs = CrossSection(eKin);
+
+  if (rj <0.0 || xs < 0.0) {  
+    std::ostringstream errOs;
+    G4cout<<"WARNING:  NEGATIVE VALUES "<<G4endl;     
+    errOs << "Rj=" << rj <<"  xsec("
+          <<eKin/MeV<<" MeV)= "<< xs <<"  A= "<<GetA()<<"  Z= "<<GetZ()
+          <<G4endl;
+    throw G4HadronicException(__FILE__, __LINE__, errOs.str());
+  }
+
+  G4double Probability = fact * GetReducedMass() * rj * xs * eKin * P * (N-1) 
+    * g4pow->powN(g1*E1/(g0*E0),N-2) * g1 / (E0*g0*g0);
   
-  G4double E0 = std::max(0.0,U - A0);
-  if (E0 == 0.0) return 0.0;
-  G4double E1 = std::max(0.0,U - eKin - GetBindingEnergy() - A1);
-  if (E1 == 0.0) return 0.0;
-
-
+  /*
   G4double Probability = 1.0/pi2*2.0/(hbarc*hbarc*hbarc) * GetReducedMass()
     * GetRj(aFragment.GetNumberOfParticles(), aFragment.GetNumberOfCharged())  
-    * eKin*CrossSection(eKin)*millibarn * P*(N-1.0) * std::pow(g1*E1/(g0*E0),N-2.0)/E0 * g1/(g0*g0);
-
-  if (GetRj(aFragment.GetNumberOfParticles(), aFragment.GetNumberOfCharged())<0.0 
-      || CrossSection(eKin) <0) {  
-    std::ostringstream errOs;
-    G4cout<<"WARNING:  NEGATIVE VALUES "<<G4endl;     
-    errOs << "Rj=" << GetRj(aFragment.GetNumberOfParticles(), aFragment.GetNumberOfCharged())
-          <<G4endl;
-    errOs <<"  xsec("<<eKin<<" MeV) ="<<CrossSection(eKin)<<G4endl;
-    errOs <<"  A="<<GetA()<<"  Z="<<GetZ()<<G4endl;
-    throw G4HadronicException(__FILE__, __LINE__, errOs.str());
-  }
+    * eKin*CrossSection(eKin)*millibarn * P*(N-1.0) * 
+   std::pow(g1*E1/(g0*E0),N-2.0)/E0 * g1/(g0*g0);
+  */
 
   return Probability;

trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/src/G4PreCompoundParameters.cc

@@ -24 +24 @@
 // ********************************************************************
 //
-//
-// $Id: G4PreCompoundParameters.cc,v 1.3 2006/06/29 20:59:29 gunter Exp $
-// GEANT4 tag $Name: geant4-09-04-beta-01 $
+// $Id: G4PreCompoundParameters.cc,v 1.4 2010/08/18 14:07:24 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03-ref-09 $
 //
 // by V. Lara
+//
+// 18.08.2010 V.Ivanchenko make this class as a standard singleton
+//
 
 #include "G4PreCompoundParameters.hh"
 
-
-G4PreCompoundParameters G4PreCompoundParameters::thePreCompoundParameters;
+G4PreCompoundParameters* G4PreCompoundParameters::theParameters = 0;
 
 G4PreCompoundParameters * G4PreCompoundParameters::GetAddress()
-{ return &thePreCompoundParameters; }
+{ 
+  if(0 == theParameters) {
+    static G4PreCompoundParameters par;
+    theParameters = ∥
+  }
+  return theParameters; 
+}
 
+G4PreCompoundParameters::G4PreCompoundParameters() 
+  : fLevelDensity(0.10/MeV),
+    fR0(1.5*fermi),
+    fTransitions_r0(0.6*fermi),
+    fFermiEnergy(35.0*MeV) 
+{}
+
+G4PreCompoundParameters::~G4PreCompoundParameters() 
+{}
+
+

trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/src/G4PreCompoundProton.cc

@@ -24 +24 @@
 // ********************************************************************
 //
-//
-// $Id: G4PreCompoundProton.cc,v 1.5 2010/04/09 14:06:17 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-04-beta-01 $
+// $Id: G4PreCompoundProton.cc,v 1.6 2010/08/28 15:16:55 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03-ref-09 $
 //
 // -------------------------------------------------------------------
@@ -39 +38 @@
 // Modified:  
 // 21.08.2008 J. M. Quesada added external choice of inverse cross section option
-// 21.08.2008 J. M. Quesada added external choice for superimposed Coulomb barrier 
-//                          (if useSICB=true) 
+// 21.08.2008 J. M. Quesada added external choice for superimposed Coulomb 
+//                          barrier (if useSICB=true) 
+// 20.08.2010 V.Ivanchenko added G4Pow and G4PreCompoundParameters pointers
+//                         use int Z and A and cleanup
 //
 
 #include "G4PreCompoundProton.hh"
-
-G4ReactionProduct * G4PreCompoundProton::GetReactionProduct() const
-{
-  G4ReactionProduct * theReactionProduct = 
-    new G4ReactionProduct(G4Proton::ProtonDefinition());
-  theReactionProduct->SetMomentum(GetMomentum().vect());
-  theReactionProduct->SetTotalEnergy(GetMomentum().e());
-#ifdef PRECOMPOUND_TEST
-  theReactionProduct->SetCreatorModel("G4PrecompoundModel");
-#endif
-  return theReactionProduct;
-}
-
-G4double G4PreCompoundProton::GetRj(const G4int NumberParticles, const G4int NumberCharged)
+#include "G4Proton.hh"
+
+G4PreCompoundProton::G4PreCompoundProton()
+  : G4PreCompoundNucleon(G4Proton::Proton(), &theProtonCoulombBarrier)
+{}
+
+G4PreCompoundProton::~G4PreCompoundProton()
+{}
+
+G4double G4PreCompoundProton::GetRj(G4int nParticles, G4int nCharged)
 {
   G4double rj = 0.0;
-  if(NumberParticles > 0) rj = static_cast<G4double>(NumberCharged)/static_cast<G4double>(NumberParticles);
+  if(nParticles > 0) { 
+    rj = static_cast<G4double>(nCharged)/static_cast<G4double>(nParticles);
+  }
   return rj;
 }
@@ -71 +70 @@
 //OPT=3 Kalbach's parameterization 
 // 
-G4double G4PreCompoundProton::CrossSection(const  G4double K)
-{
-  //G4cout<<" In G4PreCompoundProton OPTxs="<<OPTxs<<G4endl;
-  //G4cout<<" In G4PreCompoundProton useSICB="<<useSICB<<G4endl;
-
-  ResidualA=GetRestA();
-  ResidualZ=GetRestZ(); 
-  theA=GetA();
-  theZ=GetZ();
-  ResidualAthrd=std::pow(ResidualA,0.33333);
-  FragmentA=GetA()+GetRestA();
-  FragmentAthrd=std::pow(FragmentA,0.33333);
-
-  if (OPTxs==0) return GetOpt0(K);
-  else if( OPTxs==1) return GetOpt1(K);
-  else if( OPTxs==2|| OPTxs==4) return GetOpt2(K);
-  else if (OPTxs==3)  return GetOpt3(K);
+G4double G4PreCompoundProton::CrossSection(G4double K)
+{
+  ResidualA = GetRestA();
+  ResidualZ = GetRestZ(); 
+  theA = GetA();
+  theZ = GetZ();
+  ResidualAthrd = ResidualA13();
+  FragmentA = theA + ResidualA;
+  FragmentAthrd = g4pow->Z13(FragmentA);
+
+  if (OPTxs==0) { return GetOpt0(K); }
+  else if( OPTxs==1) { return GetOpt1(K); }
+  else if( OPTxs==2|| OPTxs==4) { return GetOpt2(K); }
+  else if (OPTxs==3)  { return GetOpt3(K); }
   else{
     std::ostringstream errOs;
@@ -96 +92 @@
 }
 
-// *********************** OPT=0 : Dostrovski's cross section  *****************************
-
-G4double G4PreCompoundProton::GetOpt0(const  G4double K)
-{
-  const G4double r0 = G4PreCompoundParameters::GetAddress()->Getr0();
-  // cross section is now given in mb (r0 is in mm) for the sake of consistency
-  //with the rest of the options
-  return 1.e+25*pi*(r0*ResidualAthrd)*(r0*ResidualAthrd)*GetAlpha()*(1.+GetBeta()/K);
-}
-//
-//------------
-//
 G4double G4PreCompoundProton::GetAlpha()
 {
-  G4double aZ = static_cast<G4double>(GetRestZ());
+  G4int aZ = ResidualZ;
   G4double C = 0.0;
   if (aZ >= 70) 
@@ -122 +106 @@
   return 1.0 + C;
 }
-//
-//-------------------
-//  
+
 G4double G4PreCompoundProton::GetBeta() 
 {
   return -GetCoulombBarrier();
 }
-//
-  
-//********************* OPT=1 : Chatterjee's cross section ************************ 
+  
+//********************* OPT=1 : Chatterjee's cross section *********************
 //(fitting to cross section from Bechetti & Greenles OM potential)
 
-G4double G4PreCompoundProton::GetOpt1(const  G4double K)
+G4double G4PreCompoundProton::GetOpt1(G4double K)
 {
   G4double Kc=K; 
 
   // JMQ  xsec is set constat above limit of validity
-  if (K>50)  Kc=50;
+  if (K > 50*MeV) { Kc = 50*MeV; }
 
   G4double landa, landa0, landa1, mu, mu0, mu1,nu, nu0, nu1, nu2,xs;
@@ -159 +140 @@
   p = p0 + p1/Ec + p2/(Ec*Ec);
   landa = landa0*ResidualA + landa1;
-  mu = mu0*std::pow(ResidualA,mu1);
-  nu = std::pow(ResidualA,mu1)*(nu0 + nu1*Ec + nu2*(Ec*Ec));
+
+  G4double resmu1 = g4pow->powZ(ResidualA,mu1); 
+  mu = mu0*resmu1;
+  nu = resmu1*(nu0 + nu1*Ec + nu2*(Ec*Ec));
   q = landa - nu/(Ec*Ec) - 2*p*Ec;
   r = mu + 2*nu/Ec + p*(Ec*Ec);
@@ -170 +153 @@
 
   return xs; 
-
-}
-
-//************* OPT=2 : Welisch's proton reaction cross section ************************ 
-
-G4double G4PreCompoundProton::GetOpt2(const  G4double K)
-{
-
-  G4double rnpro,rnneu,eekin,ekin,ff1,ff2,ff3,r0,fac,fac1,fac2,b0,xine_th(0);
+}
+
+//************* OPT=2 : Welisch's proton reaction cross section ***************
+
+G4double G4PreCompoundProton::GetOpt2(G4double K)
+{
+
+  G4double eekin,ekin,ff1,ff2,ff3,r0,fac,fac1,fac2,b0,xine_th(0);
  
-  //This is redundant when the Coulomb  barrier is overimposed to all cross sections 
-  //It should be kept when Coulomb barrier only imposed at OPTxs=2
-
-  if(!useSICB && K<=theCoulombBarrier) return xine_th=0.0;
+  // This is redundant when the Coulomb  barrier is overimposed to all 
+  // cross sections 
+  // It should be kept when Coulomb barrier only imposed at OPTxs=2
+
+  if(!useSICB && K<=theCoulombBarrier) { return 0.0; }
 
   eekin=K;
-  rnpro=ResidualZ;
-  rnneu=ResidualA-ResidualZ;
+  G4int rnneu=ResidualA-ResidualZ;
   ekin=eekin/1000;
   r0=1.36*1.e-15;
@@ -194 +176 @@
   fac1=b0*(1.-1./ResidualAthrd);
   fac2=1.;
-  if(rnneu > 1.5) fac2=std::log(rnneu);
+  if(rnneu > 1.5) { fac2 = g4pow->logZ(rnneu); }
   xine_th= 1.e+31*fac*fac2*(1.+ResidualAthrd-fac1);
   xine_th=(1.-0.15*std::exp(-ekin))*xine_th/(1.00-0.0007*ResidualA);    
-  ff1=0.70-0.0020*ResidualA ;
+  ff1=0.70-0.0020*ResidualA;
   ff2=1.00+1/ResidualA;
   ff3=0.8+18/ResidualA-0.002*ResidualA;
@@ -215 +197 @@
   }
   return xine_th;
-            
-}
-
+}
 
 // *********** OPT=3 : Kalbach's cross sections (from PRECO code)*************
@@ -247 +227 @@
   G4double ec,ecsq,xnulam,etest(0.),ra(0.),a,w,c,signor(1.),signor2,sig; 
   G4double b,ecut,cut,ecut2,geom,elab;
-  
-  
+    
   G4double      flow = 1.e-18;
   G4double       spill= 1.e+18; 
-  
-  
-  
-  if (ResidualA <= 60.)  signor = 0.92;
-  else if (ResidualA < 100.) signor = 0.8 + ResidualA*0.002;
-  
+   
+  if (ResidualA <= 60.)  { signor = 0.92; }
+  else if (ResidualA < 100.) { signor = 0.8 + ResidualA*0.002; }
   
   ec = 1.44 * theZ * ResidualZ / (1.5*ResidualAthrd+ra);
@@ -262 +238 @@
   p = p0 + p1/ec + p2/ecsq;
   landa = landa0*ResidualA + landa1;
-  a = std::pow(ResidualA,mu1);
+  a = g4pow->powZ(ResidualA,mu1);
   mu = mu0 * a;
   nu = a* (nu0+nu1*ec+nu2*ecsq);
@@ -270 +246 @@
   
   xnulam = nu / landa;
-  if (xnulam > spill) xnulam=0.;
-  if (xnulam >= flow) etest =std::sqrt(xnulam) + 7.;
+  if (xnulam > spill) { xnulam=0.; }
+  if (xnulam >= flow) { etest =std::sqrt(xnulam) + 7.; }
   
   a = -2.*p*ec + landa - nu/ecsq;
@@ -277 +253 @@
   ecut = 0.;
   cut = a*a - 4.*p*b;
-  if (cut > 0.) ecut = std::sqrt(cut);
+  if (cut > 0.) { ecut = std::sqrt(cut); }
   ecut = (ecut-a) / (p+p);
   ecut2 = ecut;
-//JMQ 290310 for avoiding unphysical increase below minimum (at ecut)
-//ecut<0 means that there is no cut with energy axis, i.e. xs is set to 0 bellow minimum
-//  if (cut < 0.) ecut2 = ecut - 2.;
-  if (cut < 0.) ecut2 = ecut;
+  //JMQ 290310 for avoiding unphysical increase below minimum (at ecut)
+  // ecut<0 means that there is no cut with energy axis, i.e. xs is set 
+  // to 0 bellow minimum
+  //  if (cut < 0.) ecut2 = ecut - 2.;
+  if (cut < 0.) { ecut2 = ecut; }
   elab = K * FragmentA / ResidualA;
   sig = 0.;
   if (elab <= ec) { //start for E<Ec 
-    if (elab > ecut2)  sig = (p*elab*elab+a*elab+b) * signor;
+    if (elab > ecut2) { sig = (p*elab*elab+a*elab+b) * signor; }
     
     signor2 = (ec-elab-c) / w;
@@ -308 +285 @@
     
   }   //end for E>Ec
-
   return sig;
 }
-
-//   ************************** end of cross sections ******************************* 
-
-
-

trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/src/G4PreCompoundTransitions.cc

@@ -24 +24 @@
 // ********************************************************************
 //
-// $Id: G4PreCompoundTransitions.cc,v 1.22 2009/11/21 18:03:13 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-04-beta-01 $
+// $Id: G4PreCompoundTransitions.cc,v 1.27 2010/10/20 00:47:46 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03-ref-09 $
 //
 // -------------------------------------------------------------------
@@ -37 +37 @@
 //
 // Modified:  
-// 16.02.2008 J. M. Quesada fixed bugs 
-// 06.09.2008 J. M. Quesada added external choices for:
+// 16.02.2008 J.M.Quesada fixed bugs 
+// 06.09.2008 J.M.Quesada added external choices for:
 //                      - "never go back"  hipothesis (useNGB=true) 
 //                      -  CEM transition probabilities (useCEMtr=true)
-
-// 30.10.09 J.M.Quesada: CEM transition probabilities have been renormalized 
+// 30.10.2009 J.M.Quesada: CEM transition probabilities have been renormalized 
 //                       (IAEA benchmark)
-//
+// 20.08.2010 V.Ivanchenko move constructor and destructor to the source and 
+//                         optimise the code
+//
+
 #include "G4PreCompoundTransitions.hh"
 #include "G4HadronicException.hh"
-
-const G4PreCompoundTransitions & G4PreCompoundTransitions::
-operator=(const G4PreCompoundTransitions &)
+#include "G4PreCompoundParameters.hh"
+#include "G4Proton.hh"
+#include "Randomize.hh"
+#include "G4Pow.hh"
+
+G4PreCompoundTransitions::G4PreCompoundTransitions() 
 {
-  throw G4HadronicException(__FILE__, __LINE__, "G4PreCompoundTransitions::operator= meant to not be accessable");
-  return *this;
+  proton = G4Proton::Proton();
+  FermiEnergy = G4PreCompoundParameters::GetAddress()->GetFermiEnergy();
+  r0 = G4PreCompoundParameters::GetAddress()->GetTransitionsr0();
+  aLDP = G4PreCompoundParameters::GetAddress()->GetLevelDensity();
+  g4pow = G4Pow::GetInstance();
 }
 
-
-G4bool G4PreCompoundTransitions::operator==(const G4PreCompoundTransitions &) const
-{
-  return false;
-}
-
-G4bool G4PreCompoundTransitions::operator!=(const G4PreCompoundTransitions &) const
-{
-  return true;
-}
-
-
+G4PreCompoundTransitions::~G4PreCompoundTransitions() 
+{}
+
+// Calculates transition probabilities with 
+// DeltaN = +2 (Trans1) -2 (Trans2) and 0 (Trans3)
 G4double G4PreCompoundTransitions::
 CalculateProbability(const G4Fragment & aFragment)
 {
-  //G4cout<<"In G4PreCompoundTransitions.cc  useNGB="<<useNGB<<G4endl;
-  //G4cout<<"In G4PreCompoundTransitions.cc  useCEMtr="<<useCEMtr<<G4endl;
-
-  // Fermi energy
-  const G4double FermiEnergy = G4PreCompoundParameters::GetAddress()->GetFermiEnergy();
+  // Number of holes
+  G4int H = aFragment.GetNumberOfHoles();
+  // Number of Particles 
+  G4int P = aFragment.GetNumberOfParticles();
+  // Number of Excitons 
+  G4int N = P+H;
+  // Nucleus 
+  G4int A = aFragment.GetA_asInt();
+  G4int Z = aFragment.GetZ_asInt();
+  G4double U = aFragment.GetExcitationEnergy();
+
+  //G4cout << aFragment << G4endl;
   
-  // Nuclear radius
-  const G4double r0 = G4PreCompoundParameters::GetAddress()->GetTransitionsr0();
-   
-  // In order to calculate the level density parameter
-  // G4EvaporationLevelDensityParameter theLDP;
-
-  // Number of holes
-  G4double H = aFragment.GetNumberOfHoles();
-  // Number of Particles 
-  G4double P = aFragment.GetNumberOfParticles();
-  // Number of Excitons 
-  G4double N = P+H;
-  // Nucleus 
-  G4double A = aFragment.GetA();
-  G4double Z = static_cast<G4double>(aFragment.GetZ());
-  G4double U = aFragment.GetExcitationEnergy();
-  
-  if(U<10*eV) return 0.0;
+  if(U < 10*eV) { return 0.0; }
   
   //J. M. Quesada (Feb. 08) new physics
-  // OPT=1 Transitions are calculated according to Gudima's paper (original in G4PreCompound from VL) 
+  // OPT=1 Transitions are calculated according to Gudima's paper 
+  //       (original in G4PreCompound from VL) 
   // OPT=2 Transitions are calculated according to Gupta's formulae
   //
-  
-  
-  
   if (useCEMtr){
 
-    
     // Relative Energy (T_{rel})
-    G4double RelativeEnergy = (8.0/5.0)*FermiEnergy + U/N;
+    G4double RelativeEnergy = 1.6*FermiEnergy + U/G4double(N);
     
     // Sample kind of nucleon-projectile 
     G4bool ChargedNucleon(false);
     G4double chtest = 0.5;
-    if (P > 0) chtest = aFragment.GetNumberOfCharged()/P;
-    if (G4UniformRand() < chtest) ChargedNucleon = true;
+    if (P > 0) { 
+      chtest = G4double(aFragment.GetNumberOfCharged())/G4double(P); 
+    }
+    if (G4UniformRand() < chtest) { ChargedNucleon = true; }
     
     // Relative Velocity: 
     // <V_{rel}>^2
     G4double RelativeVelocitySqr(0.0);
-    if (ChargedNucleon) RelativeVelocitySqr = 2.0*RelativeEnergy/proton_mass_c2;
-    else RelativeVelocitySqr = 2.0*RelativeEnergy/neutron_mass_c2;
+    if (ChargedNucleon) { 
+      RelativeVelocitySqr = 2.0*RelativeEnergy/CLHEP::proton_mass_c2; 
+    } else { 
+      RelativeVelocitySqr = 2.0*RelativeEnergy/CLHEP::neutron_mass_c2; 
+    }
     
     // <V_{rel}>
@@ -124 +117 @@
     
     // Proton-Proton Cross Section
-    G4double ppXSection = (10.63/RelativeVelocitySqr - 29.92/RelativeVelocity + 42.9)*millibarn;
+    G4double ppXSection = 
+      (10.63/RelativeVelocitySqr - 29.92/RelativeVelocity + 42.9)
+      * CLHEP::millibarn;
     // Proton-Neutron Cross Section
-    G4double npXSection = (34.10/RelativeVelocitySqr - 82.20/RelativeVelocity + 82.2)*millibarn;
+    G4double npXSection = 
+      (34.10/RelativeVelocitySqr - 82.20/RelativeVelocity + 82.2)
+      * CLHEP::millibarn;
     
     // Averaged Cross Section: \sigma(V_{rel})
@@ -134 +131 @@
       {
         //JMQ: small bug fixed
-        //      AveragedXSection = ((Z-1.0) * ppXSection + (A-Z-1.0) * npXSection) / (A-1.0);
-        AveragedXSection = ((Z-1.0) * ppXSection + (A-Z) * npXSection) / (A-1.0);
+        //AveragedXSection=((Z-1.0) * ppXSection + (A-Z-1.0) * npXSection)/(A-1.0);
+        AveragedXSection = ((Z-1)*ppXSection + (A-Z)*npXSection)/G4double(A-1);
       }
     else 
       {
-        AveragedXSection = ((A-Z-1.0) * ppXSection + Z * npXSection) / (A-1.0);
+        AveragedXSection = ((A-Z-1)*ppXSection + Z*npXSection)/G4double(A-1);
+        //AveragedXSection = ((A-Z-1)*npXSection + Z*ppXSection)/G4double(A-1);
       }
     
@@ -146 +144 @@
     
     // This factor is introduced to take into account the Pauli principle
-    G4double PauliFactor = 1.0 - (7.0/5.0)*FermiRelRatio;
-    if (FermiRelRatio > 0.5) PauliFactor += (2.0/5.0)*FermiRelRatio*std::pow(2.0 - (1.0/FermiRelRatio), 5.0/2.0);
-    
+    G4double PauliFactor = 1.0 - 1.4*FermiRelRatio;
+    if (FermiRelRatio > 0.5) {
+      G4double x = 2.0 - 1.0/FermiRelRatio;
+      PauliFactor += 0.4*FermiRelRatio*x*x*std::sqrt(x);
+      //PauliFactor += 
+      //(2.0/5.0)*FermiRelRatio*std::pow(2.0 - (1.0/FermiRelRatio), 5.0/2.0);
+    }
     // Interaction volume 
-    //  G4double Vint = (4.0/3.0)*pi*std::pow(2.0*r0 + hbarc/(proton_mass_c2*RelativeVelocity) , 3.0);
-    G4double xx=2.0*r0 + hbarc/(proton_mass_c2*RelativeVelocity);
-    G4double Vint = (4.0/3.0)*pi*xx*xx*xx;
+    //  G4double Vint = (4.0/3.0)
+    //*pi*std::pow(2.0*r0 + hbarc/(proton_mass_c2*RelativeVelocity) , 3.0);
+    G4double xx = 2.0*r0 + hbarc/(CLHEP::proton_mass_c2*RelativeVelocity);
+    //    G4double Vint = (4.0/3.0)*CLHEP::pi*xx*xx*xx;
+    G4double Vint = CLHEP::pi*xx*xx*xx/0.75;
     
     // Transition probability for \Delta n = +2
     
-    TransitionProb1 = AveragedXSection*PauliFactor*std::sqrt(2.0*RelativeEnergy/proton_mass_c2)/Vint;
-
-//JMQ 281009  phenomenological factor in order to increase equilibrium contribution
-//   G4double factor=5.0;
-//   TransitionProb1 *= factor;
-//
-    if (TransitionProb1 < 0.0) TransitionProb1 = 0.0; 
-    
-    G4double a = G4PreCompoundParameters::GetAddress()->GetLevelDensity();
+    TransitionProb1 = AveragedXSection*PauliFactor
+      *std::sqrt(2.0*RelativeEnergy/CLHEP::proton_mass_c2)/Vint;
+
+    //JMQ 281009  phenomenological factor in order to increase 
+    //   equilibrium contribution
+    //   G4double factor=5.0;
+    //   TransitionProb1 *= factor;
+    //
+    if (TransitionProb1 < 0.0) { TransitionProb1 = 0.0; } 
+    
     // GE = g*E where E is Excitation Energy
-    G4double GE = (6.0/pi2)*a*A*U;
-    
-    G4double Fph = ((P*P+H*H+P-H)/4.0 - H/2.0);
-    
-    //G4bool NeverGoBack(false);
-    G4bool NeverGoBack;
-    if(useNGB)  NeverGoBack=true;
-    else NeverGoBack=false;
-    
+    G4double GE = (6.0/pi2)*aLDP*A*U;
+    
+    //G4double Fph = ((P*P+H*H+P-H)/4.0 - H/2.0);
+    G4double Fph = G4double(P*P+H*H+P-3*H)/4.0;
+    
+    G4bool NeverGoBack(false);
+    if(useNGB) { NeverGoBack=true; }
     
     //JMQ/AH  bug fixed: if (U-Fph < 0.0) NeverGoBack = true;
-    if (GE-Fph < 0.0) NeverGoBack = true;
+    if (GE-Fph < 0.0) { NeverGoBack = true; }
     
     // F(p+1,h+1)
     G4double Fph1 = Fph + N/2.0;
     
-    G4double ProbFactor = std::pow((GE-Fph)/(GE-Fph1),N+1.0);
-    
+    G4double ProbFactor = g4pow->powN((GE-Fph)/(GE-Fph1),N+1);
     
     if (NeverGoBack)
       {
-      TransitionProb2 = 0.0;
-      TransitionProb3 = 0.0;
+        TransitionProb2 = 0.0;
+        TransitionProb3 = 0.0;
       }
     else 
       {
         // Transition probability for \Delta n = -2 (at F(p,h) = 0)
-        TransitionProb2 = TransitionProb1 * ProbFactor * (P*H*(N+1.0)*(N-2.0))/((GE-Fph)*(GE-Fph));
-        if (TransitionProb2 < 0.0) TransitionProb2 = 0.0; 
+        TransitionProb2 = 
+          TransitionProb1 * ProbFactor * (P*H*(N+1)*(N-2))/((GE-Fph)*(GE-Fph));
+        if (TransitionProb2 < 0.0) { TransitionProb2 = 0.0; } 
         
         // Transition probability for \Delta n = 0 (at F(p,h) = 0)
-        TransitionProb3 = TransitionProb1* ((N+1.0)/N) * ProbFactor  * (P*(P-1.0) + 4.0*P*H + H*(H-1.0))/(GE-Fph);
-        if (TransitionProb3 < 0.0) TransitionProb3 = 0.0; 
-      }
-    
-    //  G4cout<<"U = "<<U<<G4endl;
-    //  G4cout<<"N="<<N<<"  P="<<P<<"  H="<<H<<G4endl;
-    //  G4cout<<"l+ ="<<TransitionProb1<<"  l- ="<< TransitionProb2<<"  l0 ="<< TransitionProb3<<G4endl; 
-    return TransitionProb1 + TransitionProb2 + TransitionProb3;
-  }
-  
-  else  {
-    //JMQ: Transition probabilities from Gupta's work
-    
-    G4double a = G4PreCompoundParameters::GetAddress()->GetLevelDensity();
+        TransitionProb3 = TransitionProb1*(N+1)* ProbFactor  
+          * (P*(P-1) + 4.0*P*H + H*(H-1))/(N*(GE-Fph));
+        if (TransitionProb3 < 0.0) { TransitionProb3 = 0.0; }
+      }
+  } else {
+    //JMQ: Transition probabilities from Gupta's work    
     // GE = g*E where E is Excitation Energy
-    G4double GE = (6.0/pi2)*a*A*U;
+    G4double GE = (6.0/pi2)*aLDP*A*U;
  
     G4double Kmfp=2.;
         
-    TransitionProb1=1./Kmfp*3./8.*1./c_light*1.0e-9*(1.4e+21*U-2./(N+1)*6.0e+18*U*U);
-    if (TransitionProb1 < 0.0) TransitionProb1 = 0.0;
-    
-    if (useNGB){
-      TransitionProb2=0.;
-      TransitionProb3=0.;
-    }
-    else{        
-      if (N<=1) TransitionProb2=0. ;    
-      else  TransitionProb2=1./Kmfp*3./8.*1./c_light*1.0e-9*(N-1.)*(N-2.)*P*H/(GE*GE)*(1.4e+21*U - 2./(N-1)*6.0e+18*U*U);      
+    //TransitionProb1=1./Kmfp*3./8.*1./c_light*1.0e-9*(1.4e+21*U-2./(N+1)*6.0e+18*U*U);
+    TransitionProb1 = 3.0e-9*(1.4e+21*U - 1.2e+19*U*U/G4double(N+1))
+      /(8*Kmfp*CLHEP::c_light);
+    if (TransitionProb1 < 0.0) { TransitionProb1 = 0.0; }
+
+    TransitionProb2=0.;
+    TransitionProb3=0.;
+    
+    if (!useNGB && N > 1) {
+      // TransitionProb2=1./Kmfp*3./8.*1./c_light*1.0e-9*(N-1.)*(N-2.)*P*H/(GE*GE)*(1.4e+21*U - 2./(N-1)*6.0e+18*U*U);      
+      TransitionProb2 = 
+        3.0e-9*(N-2)*P*H*(1.4e+21*U*(N-1) - 1.2e+19*U*U)/(8*Kmfp*c_light*GE*GE);      
       if (TransitionProb2 < 0.0) TransitionProb2 = 0.0; 
-      TransitionProb3=0.;
-    }
-    
-      //  G4cout<<"U = "<<U<<G4endl;
-    //  G4cout<<"N="<<N<<"  P="<<P<<"  H="<<H<<G4endl;
-    //  G4cout<<"l+ ="<<TransitionProb1<<"  l- ="<< TransitionProb2<<"  l0 ="<< TransitionProb3<<G4endl; 
-    return TransitionProb1 + TransitionProb2 + TransitionProb3;
+    }
   }
+  //  G4cout<<"U = "<<U<<G4endl;
+  //  G4cout<<"N="<<N<<"  P="<<P<<"  H="<<H<<G4endl;
+  //  G4cout<<"l+ ="<<TransitionProb1<<"  l- ="<< TransitionProb2
+  //   <<"  l0 ="<< TransitionProb3<<G4endl; 
+  return TransitionProb1 + TransitionProb2 + TransitionProb3;
 }
 
-
-G4Fragment G4PreCompoundTransitions::PerformTransition(const G4Fragment & aFragment)
+void G4PreCompoundTransitions::PerformTransition(G4Fragment & result)
 {
-  G4Fragment result(aFragment);
-  G4double ChosenTransition = G4UniformRand()*(TransitionProb1 + TransitionProb2 + TransitionProb3);
+  G4double ChosenTransition = 
+    G4UniformRand()*(TransitionProb1 + TransitionProb2 + TransitionProb3);
   G4int deltaN = 0;
-  G4int Nexcitons = result.GetNumberOfExcitons();
+  //  G4int Nexcitons = result.GetNumberOfExcitons();
+  G4int Npart     = result.GetNumberOfParticles();
+  G4int Ncharged  = result.GetNumberOfCharged();
+  G4int Nholes    = result.GetNumberOfHoles();
   if (ChosenTransition <= TransitionProb1) 
     {
@@ -255 +253 @@
     }
 
-  // AH/JMQ: Randomly decrease the number of charges if deltaN is -2 and in proportion 
-  // to the number charges w.r.t. number of particles,  PROVIDED that there are charged particles
-  if(deltaN < 0 && G4UniformRand() <= 
-     static_cast<G4double>(result.GetNumberOfCharged())/static_cast<G4double>(result.GetNumberOfParticles()) 
-     && (result.GetNumberOfCharged() >= 1)) {
-    result.SetNumberOfCharged(result.GetNumberOfCharged()+deltaN/2); // deltaN is negative!
-  }
+  // AH/JMQ: Randomly decrease the number of charges if deltaN is -2 and  
+  // in proportion to the number charges w.r.t. number of particles,  
+  // PROVIDED that there are charged particles
+  deltaN /= 2;
+
+  //G4cout << "deltaN= " << deltaN << G4endl;
 
   // JMQ the following lines have to be before SetNumberOfCharged, otherwise the check on 
   // number of charged vs. number of particles fails
-  result.SetNumberOfParticles(result.GetNumberOfParticles()+deltaN/2);
-  result.SetNumberOfHoles(result.GetNumberOfHoles()+deltaN/2); 
-
-  // With weight Z/A, number of charged particles is increased with +1
-  if ( ( deltaN > 0 ) &&
-      (G4UniformRand() <= static_cast<G4double>(result.GetZ()-result.GetNumberOfCharged())/
-       std::max(static_cast<G4double>(result.GetA()-Nexcitons),1.)))
-    {
-      result.SetNumberOfCharged(result.GetNumberOfCharged()+deltaN/2);
-    }
+  result.SetNumberOfParticles(Npart+deltaN);
+  result.SetNumberOfHoles(Nholes+deltaN); 
+
+  if(deltaN < 0) {
+    if( Ncharged >= 1 && G4int(Npart*G4UniformRand()) <= Ncharged) 
+      { 
+        result.SetNumberOfCharged(Ncharged+deltaN); // deltaN is negative!
+      }
+
+  } else if ( deltaN > 0 ) {
+    // With weight Z/A, number of charged particles is increased with +1
+    G4int A = result.GetA_asInt();
+    G4int Z = result.GetZ_asInt();
+    if( G4int(std::max(1, A - Npart)*G4UniformRand()) <= Z) 
+      {
+        result.SetNumberOfCharged(Ncharged+deltaN);
+      }
+  }
   
   // Number of charged can not be greater that number of particles
-  if ( result.GetNumberOfParticles() < result.GetNumberOfCharged() ) 
+  if ( Npart < Ncharged ) 
     {
-      result.SetNumberOfCharged(result.GetNumberOfParticles());
-    }
-  
-  return result;
+      result.SetNumberOfCharged(Npart);
+    }
+  //G4cout << "### After transition" << G4endl;
+  //G4cout << result << G4endl;
 }
 

trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/src/G4PreCompoundTriton.cc

@@ -24 +24 @@
 // ********************************************************************
 //
-//
-// $Id: G4PreCompoundTriton.cc,v 1.6 2010/04/09 14:06:17 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-04-beta-01 $
+// $Id: G4PreCompoundTriton.cc,v 1.7 2010/08/28 15:16:55 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03-ref-09 $
 //
 // -------------------------------------------------------------------
@@ -39 +38 @@
 // Modified:  
 // 21.08.2008 J. M. Quesada add choice of options  
+// 20.08.2010 V.Ivanchenko added G4Pow and G4PreCompoundParameters pointers
+//                         use int Z and A and cleanup
 //
  
 #include "G4PreCompoundTriton.hh"
-
-
-G4ReactionProduct * G4PreCompoundTriton::GetReactionProduct() const
-{
-  G4ReactionProduct * theReactionProduct =
-    new G4ReactionProduct(G4Triton::TritonDefinition());
-  theReactionProduct->SetMomentum(GetMomentum().vect());
-  theReactionProduct->SetTotalEnergy(GetMomentum().e());
-#ifdef PRECOMPOUND_TEST
-  theReactionProduct->SetCreatorModel("G4PrecompoundModel");
-#endif
-  return theReactionProduct;
-}   
-
-G4double G4PreCompoundTriton::FactorialFactor(const G4double N, const G4double P)
-{
-  return 
-      (N-3.0)*(P-2.0)*(
-                       (((N-2.0)*(P-1.0))/2.0) *(
-                                                 (((N-1.0)*P)/3.0) 
-                                                 )
-                       );
+#include "G4Triton.hh"
+
+G4PreCompoundTriton::G4PreCompoundTriton()
+  : G4PreCompoundIon(G4Triton::Triton(), &theTritonCoulombBarrier)
+{}
+
+G4PreCompoundTriton::~G4PreCompoundTriton()
+{}
+
+G4double G4PreCompoundTriton::FactorialFactor(G4int N, const G4int P)
+{
+  return G4double((N-3)*(P-2)*(N-2)*(P-1)*(N-1)*P)/6.0; 
 }
   
-G4double G4PreCompoundTriton::CoalescenceFactor(const G4double A)
-{
-  return 243.0/(A*A);
+G4double G4PreCompoundTriton::CoalescenceFactor(G4int A)
+{
+  return 243.0/G4double(A*A);
 }    
 
-G4double G4PreCompoundTriton::GetRj(const G4int NumberParticles, const G4int NumberCharged)
+G4double G4PreCompoundTriton::GetRj(G4int nParticles, G4int nCharged)
 {
   G4double rj = 0.0;
-  G4double denominator = NumberParticles*(NumberParticles-1)*(NumberParticles-2);
-  if(NumberCharged >= 1 && (NumberParticles-NumberCharged) >= 2) {
-    rj = 3.0*static_cast<G4double>(NumberCharged*(NumberParticles-NumberCharged)*(NumberParticles-NumberCharged-1))
-      /static_cast<G4double>(denominator); 
+  if(nCharged >= 1 && (nParticles-nCharged) >= 2) {
+    G4double denominator = 
+      G4double(nParticles*(nParticles-1)*(nParticles-2));
+    rj = G4double(3*nCharged*(nParticles-nCharged)*(nParticles-nCharged-1))
+      /denominator; 
   }
   return rj;
 }
 
-////////////////////////////////////////////////////////////////////////////////////
+//////////////////////////////////////////////////////////////////////////////////
 //J. M. Quesada (Dec 2007-June 2008): New inverse reaction cross sections 
 //OPT=0 Dostrovski's parameterization
@@ -88 +80 @@
 //OPT=3,4 Kalbach's parameterization 
 // 
-G4double G4PreCompoundTriton::CrossSection(const  G4double K)
-{
-  ResidualA=GetRestA();
-  ResidualZ=GetRestZ(); 
-  theA=GetA();
-  theZ=GetZ();
-  ResidualAthrd=std::pow(ResidualA,0.33333);
-  FragmentA=GetA()+GetRestA();
-  FragmentAthrd=std::pow(FragmentA,0.33333);
-
-  if (OPTxs==0) return GetOpt0( K);
-  else if( OPTxs==1 || OPTxs==2) return GetOpt12( K);
-  else if (OPTxs==3 || OPTxs==4)  return GetOpt34( K);
+G4double G4PreCompoundTriton::CrossSection(G4double K)
+{
+  ResidualA = GetRestA();
+  ResidualZ = GetRestZ(); 
+  theA = GetA();
+  theZ = GetZ();
+  ResidualAthrd = ResidualA13();
+  FragmentA = theA + ResidualA;
+  FragmentAthrd = g4pow->Z13(FragmentA);
+
+  if (OPTxs==0) { return GetOpt0( K); }
+  else if( OPTxs==1 || OPTxs==2) { return GetOpt12( K); }
+  else if (OPTxs==3 || OPTxs==4) { return GetOpt34( K); }
   else{
     std::ostringstream errOs;
@@ -109 +101 @@
 }
 
-// *********************** OPT=0 : Dostrovski's cross section  *****************************
-
-G4double G4PreCompoundTriton::GetOpt0(const  G4double K)
-{
-  const G4double r0 = G4PreCompoundParameters::GetAddress()->Getr0();
-  // cross section is now given in mb (r0 is in mm) for the sake of consistency
-  //with the rest of the options
-  return 1.e+25*pi*(r0*ResidualAthrd)*(r0*ResidualAthrd)*GetAlpha()*(1.+GetBeta()/K);
-}
-//
-//---------
-//
 G4double G4PreCompoundTriton::GetAlpha()
 {
   G4double C = 0.0;
-  G4double aZ = GetZ() + GetRestZ();
+  G4int aZ = theZ + ResidualZ;
   if (aZ >= 70) 
     {
@@ -136 +116 @@
   return 1.0 + C/3.0;
 }
-//
-//-------------
-//
-G4double G4PreCompoundTriton::GetBeta() 
-{
-  return -GetCoulombBarrier();
-}
-//
-//********************* OPT=1,2 : Chatterjee's cross section ************************ 
+
+//
+//********************* OPT=1,2 : Chatterjee's cross section *****************
 //(fitting to cross section from Bechetti & Greenles OM potential)
 
-G4double G4PreCompoundTriton::GetOpt12(const  G4double K)
-{
-
+G4double G4PreCompoundTriton::GetOpt12(G4double K)
+{
   G4double Kc=K;
 
   // JMQ xsec is set constat above limit of validity
-  if (K>50) Kc=50;
+  if (K > 50*MeV) { Kc=50*MeV; }
 
   G4double landa ,mu ,nu ,p , Ec,q,r,ji,xs;
@@ -172 +145 @@
   p = p0 + p1/Ec + p2/(Ec*Ec);
   landa = landa0*ResidualA + landa1;
-  mu = mu0*std::pow(ResidualA,mu1);
-  nu = std::pow(ResidualA,mu1)*(nu0 + nu1*Ec + nu2*(Ec*Ec));
+
+  G4double resmu1 = g4pow->powZ(ResidualA,mu1); 
+  mu = mu0*resmu1;
+  nu = resmu1*(nu0 + nu1*Ec + nu2*(Ec*Ec));
   q = landa - nu/(Ec*Ec) - 2*p*Ec;
   r = mu + 2*nu/Ec + p*(Ec*Ec);
@@ -184 +159 @@
               
   return xs;
-
 }
 
 // *********** OPT=3,4 : Kalbach's cross sections (from PRECO code)*************
-G4double G4PreCompoundTriton::GetOpt34(const  G4double K)
+G4double G4PreCompoundTriton::GetOpt34(G4double K)
 //     ** t from o.m. of hafele, flynn et al
 {
-
   G4double landa, mu, nu, p , signor(1.),sig;
   G4double ec,ecsq,xnulam,etest(0.),a; 
   G4double b,ecut,cut,ecut2,geom,elab;
 
-
   G4double     flow = 1.e-18;
   G4double     spill= 1.e+18;
-
 
   G4double     p0 = -21.45;
@@ -220 +191 @@
   p = p0 + p1/ec + p2/ecsq;
   landa = landa0*ResidualA + landa1;
-  a = std::pow(ResidualA,mu1);
+  a = g4pow->powZ(ResidualA,mu1);
   mu = mu0 * a;
   nu = a* (nu0+nu1*ec+nu2*ecsq);  
   xnulam = nu / landa;
-  if (xnulam > spill) xnulam=0.;
-  if (xnulam >= flow) etest = 1.2 *std::sqrt(xnulam);
+  if (xnulam > spill) { xnulam=0.; }
+  if (xnulam >= flow) { etest = 1.2 *std::sqrt(xnulam); }
  
   a = -2.*p*ec + landa - nu/ecsq;
@@ -231 +202 @@
   ecut = 0.;
   cut = a*a - 4.*p*b;
-  if (cut > 0.) ecut = std::sqrt(cut);
+  if (cut > 0.) { ecut = std::sqrt(cut); }
   ecut = (ecut-a) / (p+p);
   ecut2 = ecut;
-//JMQ 290310 for avoiding unphysical increase below minimum (at ecut)
-//ecut<0 means that there is no cut with energy axis, i.e. xs is set to 0 bellow minimum
-//  if (cut < 0.) ecut2 = ecut - 2.;
-  if (cut < 0.) ecut2 = ecut;
-  elab = K * FragmentA / ResidualA;
+  //JMQ 290310 for avoiding unphysical increase below minimum (at ecut)
+  // ecut<0 means that there is no cut with energy axis, i.e. xs is set 
+  // to 0 bellow minimum
+  //  if (cut < 0.) ecut2 = ecut - 2.;
+  if (cut < 0.) { ecut2 = ecut; }
+  elab = K * FragmentA / G4double(ResidualA);
   sig = 0.;
  
   if (elab <= ec) { //start for E<Ec
-    if (elab > ecut2)  sig = (p*elab*elab+a*elab+b) * signor;
+    if (elab > ecut2) { sig = (p*elab*elab+a*elab+b) * signor; }
   }           //end for E<Ec
   else {           //start for E>Ec
     sig = (landa*elab+mu+nu/elab) * signor;
     geom = 0.;
-    if (xnulam < flow || elab < etest) return sig;
+    if (xnulam < flow || elab < etest) { return sig; }
     geom = std::sqrt(theA*K);
     geom = 1.23*ResidualAthrd + ra + 4.573/geom;
@@ -254 +226 @@
   }           //end for E>Ec
   return sig;
-
-}
-
-//   ************************** end of cross sections ******************************* 
+}

trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/src/G4VPreCompoundEmissionFactory.cc

@@ -24 +24 @@
 // ********************************************************************
 //
-
+// $Id: G4VPreCompoundEmissionFactory.cc,v 1.5 2010/08/28 15:16:55 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03-ref-09 $
+//
 // by V. Lara
 
 #include "G4VPreCompoundEmissionFactory.hh"
-#include "G4HadronicException.hh"
 
-const G4VPreCompoundEmissionFactory & 
-G4VPreCompoundEmissionFactory::operator=(const G4VPreCompoundEmissionFactory & )
+G4VPreCompoundEmissionFactory::G4VPreCompoundEmissionFactory() 
+  : fragvector(0) 
+{}
+
+G4VPreCompoundEmissionFactory::~G4VPreCompoundEmissionFactory()
 {
-  throw G4HadronicException(__FILE__, __LINE__, "G4VPreCompoundEmissionFactory::operator= meant to not be accessable.");
-  return *this;
-}
-
-G4bool 
-G4VPreCompoundEmissionFactory::operator==(const G4VPreCompoundEmissionFactory & ) const
-{
-  throw G4HadronicException(__FILE__, __LINE__, "G4VPreCompoundEmissionFactory::operator== meant to not be accessable.");
-  return false;
-}
-
-G4bool 
-G4VPreCompoundEmissionFactory::operator!=(const G4VPreCompoundEmissionFactory & ) const
-{
-  throw G4HadronicException(__FILE__, __LINE__, "G4VPreCompoundEmissionFactory::operator!= meant to not be accessable.");
-  return true;
+  if (fragvector != 0)
+    std::for_each(fragvector->begin(), fragvector->end(), 
+                  DeleteFragment());
+  delete fragvector;
 }
 
 
-
-
-G4VPreCompoundEmissionFactory::~G4VPreCompoundEmissionFactory()
-{
-  if (_fragvector != 0)
-    std::for_each(_fragvector->begin(), _fragvector->end(), 
-                    DeleteFragment());
-  delete _fragvector;
-}
-
-
-std::vector<G4VPreCompoundFragment*> * 
-G4VPreCompoundEmissionFactory::GetFragmentVector()
-{
-  // Lazy initialization
-  if (_fragvector == 0)
-    _fragvector = CreateFragmentVector();
-  return _fragvector;
-}
-

trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/src/G4VPreCompoundFragment.cc

@@ -24 +24 @@
 // ********************************************************************
 //
-// $Id: G4VPreCompoundFragment.cc,v 1.12 2009/02/10 16:01:37 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-04-beta-01 $
+// $Id: G4VPreCompoundFragment.cc,v 1.13 2010/08/28 15:16:55 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03-ref-09 $
 //
-// J. M. Quesada (August 2008).  
-// Based  on previous work by V. Lara
+// J. M. Quesada (August 2008).  Based  on previous work by V. Lara
 //
- 
+// Modified:
+// 20.08.2010 V.Ivanchenko added G4Pow and G4PreCompoundParameters pointers
+//                         use int Z and A and cleanup
+
 #include "G4VPreCompoundFragment.hh"
 #include "G4PreCompoundParameters.hh"
+#include "G4NucleiProperties.hh"
 
-G4VPreCompoundFragment::
-G4VPreCompoundFragment(const G4VPreCompoundFragment & right)
+G4VPreCompoundFragment::G4VPreCompoundFragment(
+  const G4ParticleDefinition* part, G4VCoulombBarrier* aCoulombBarrier)
+  : particle(part), theCoulombBarrierPtr(aCoulombBarrier),
+    theRestNucleusA(0),theRestNucleusZ(0),theBindingEnergy(0.0), 
+    theMaximalKineticEnergy(-MeV),theRestNucleusMass(0.0),
+    theReducedMass(0.0),theMomentum(0.,0.,0.,0.),
+    theEmissionProbability(0.0),theCoulombBarrier(0.0)
 {
-  theA = right.theA;
-  theZ = right.theZ;
-  theRestNucleusA = right.theRestNucleusA;
-  theRestNucleusZ = right.theRestNucleusZ;
-  theCoulombBarrier = right.theCoulombBarrier;
-  theCoulombBarrierPtr = right.theCoulombBarrierPtr;
-  theMaximalKineticEnergy = right.theMaximalKineticEnergy;
-  theEmissionProbability = right.theEmissionProbability;
-  theMomentum = right.theMomentum;
-  theFragmentName = right.theFragmentName;
-  theStage = right.theStage;
+  theA = particle->GetBaryonNumber();
+  theZ = G4int(particle->GetPDGCharge()/eplus + 0.1);
+  theMass = particle->GetPDGMass();
+  theParameters = G4PreCompoundParameters::GetAddress();
+  g4pow = G4Pow::GetInstance();
 }
 
-G4VPreCompoundFragment::
-G4VPreCompoundFragment(const G4double anA,
-                       const G4double aZ, 
-                       G4VCoulombBarrier* aCoulombBarrier,
-                       const G4String & aName):
-  theA(anA),theZ(aZ), 
-  theRestNucleusA(0.0),theRestNucleusZ(0.0),theCoulombBarrier(0.0),
-  theCoulombBarrierPtr(aCoulombBarrier),
-  theBindingEnergy(0.0), theMaximalKineticEnergy(-1.0),
-  theEmissionProbability(0.0), theMomentum(0.0,0.0,0.0,0.0),
-  theFragmentName(aName),theStage(0)
+G4VPreCompoundFragment::~G4VPreCompoundFragment()
 {}
-
-
-
-G4VPreCompoundFragment::~G4VPreCompoundFragment()
-{
-}
-
-
-const G4VPreCompoundFragment & G4VPreCompoundFragment::
-operator= (const G4VPreCompoundFragment & right)
-{
-  if (this != &right) {
-    theA = right.theA;
-    theZ = right.theZ;
-    theRestNucleusA = right.theRestNucleusA;
-    theRestNucleusZ = right.theRestNucleusZ;
-    theCoulombBarrier = right.theCoulombBarrier;
-    theCoulombBarrierPtr = right.theCoulombBarrierPtr;
-    theMaximalKineticEnergy = right.theMaximalKineticEnergy;
-    theEmissionProbability = right.theEmissionProbability;
-    theMomentum = right.theMomentum;
-    theFragmentName = right.theFragmentName;
-    theStage = right.theStage;
-  }
-  return *this;
-}
-
-G4int G4VPreCompoundFragment::operator==(const G4VPreCompoundFragment & right) const
-{
-  return (this == (G4VPreCompoundFragment *) &right);
-}
-
-G4int G4VPreCompoundFragment::operator!=(const G4VPreCompoundFragment & right) const
-{
-  return (this != (G4VPreCompoundFragment *) &right);
-}
-
 
 std::ostream& 
@@ -107 +62 @@
 }
 
-
 std::ostream& 
 operator << (std::ostream &out, const G4VPreCompoundFragment *theFragment)
@@ -115 +69 @@
     
   out 
-    << "PreCompound Model Emitted Fragment: A = " 
+    << "PreCompoundModel Emitted Fragment: A = " 
     << std::setprecision(3) << theFragment->theA 
     << ", Z = " << std::setprecision(3) << theFragment->theZ;
@@ -130 +84 @@
     
     out.setf(old_floatfield,std::ios::floatfield);
-    
     return out;
 }
-
 
 void G4VPreCompoundFragment::
 Initialize(const G4Fragment & aFragment)
 {
-  theRestNucleusA = aFragment.GetA() - theA;
-  theRestNucleusZ = aFragment.GetZ() - theZ;
+  theRestNucleusA = aFragment.GetA_asInt() - theA;
+  theRestNucleusZ = aFragment.GetZ_asInt() - theZ;
+  theRestNucleusA13 = g4pow->Z13(theRestNucleusA);
 
   if ((theRestNucleusA < theRestNucleusZ) ||
@@ -149 +102 @@
       return;
     }
-  
-  
+    
   // Calculate Coulomb barrier
   theCoulombBarrier = theCoulombBarrierPtr->
-    GetCoulombBarrier(static_cast<G4int>(theRestNucleusA),static_cast<G4int>(theRestNucleusZ),
+    GetCoulombBarrier(theRestNucleusA,theRestNucleusZ,
                       aFragment.GetExcitationEnergy());
 
+  // Calculate masses
+  theRestNucleusMass = 
+    G4NucleiProperties::GetNuclearMass(theRestNucleusA, theRestNucleusZ);
+  theReducedMass = theRestNucleusMass*theMass/(theRestNucleusMass + theMass);
 
   // Compute Binding Energies for fragments 
   // (needed to separate a fragment from the nucleus)
+  theBindingEnergy = theRestNucleusMass + theMass - aFragment.GetGroundStateMass();
   
-  theBindingEnergy = G4NucleiProperties::GetMassExcess(static_cast<G4int>(theA),static_cast<G4int>(theZ)) +
-    G4NucleiProperties::GetMassExcess(static_cast<G4int>(theRestNucleusA),static_cast<G4int>(theRestNucleusZ)) -
-    G4NucleiProperties::GetMassExcess(static_cast<G4int>(aFragment.GetA()),static_cast<G4int>(aFragment.GetZ()));
+  //theBindingEnergy = G4NucleiProperties::GetMassExcess(static_cast<G4int>(theA),static_cast<G4int>(theZ)) +
+  //G4NucleiProperties::GetMassExcess(static_cast<G4int>(theRestNucleusA),static_cast<G4int>(theRestNucleusZ)) -
+  //G4NucleiProperties::GetMassExcess(static_cast<G4int>(aFragment.GetA()),static_cast<G4int>(aFragment.GetZ()));
   
   // Compute Maximal Kinetic Energy which can be carried by fragments after separation
   // This is the true (assimptotic) maximal kinetic energy
-  G4double m = aFragment.GetMomentum().m();
-  G4double rm = GetRestNuclearMass();
-  G4double em = GetNuclearMass();
+  G4double m  = aFragment.GetMomentum().m();
+  G4double rm = theRestNucleusMass;
+  G4double em = theMass;
   theMaximalKineticEnergy = ((m - rm)*(m + rm) + em*em)/(2.0*m) - em;
  
-  
   return;
 }
-
-
-
-