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

Timestamp:

Apr 6, 2009, 12:30:29 PM (15 years ago)

Author:

garnier

Message:

update processes

File:

• trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/src/G4PreCompoundNucleon.cc (modified) (4 diffs)

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

@@ -24 +24 @@
 // ********************************************************************
 //
-// by V. Lara
+//
+// $Id: G4PreCompoundNucleon.cc,v 1.13 2009/02/11 18:06:00 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-02-ref-02 $
+//
+// -------------------------------------------------------------------
+//
+// GEANT4 Class file
+//
+//
+// File name:     G4PreCompoundNucleon
+//
+// Author:         V.Lara
+//
+// Modified:  
+// 10.02.2009 J. M. Quesada cleanup
+//
 
 #include "G4PreCompoundNucleon.hh"
 #include "G4PreCompoundParameters.hh"
 
+G4bool G4PreCompoundNucleon::IsItPossible(const G4Fragment& aFragment) 
+{
+  G4int pplus = aFragment.GetNumberOfCharged();   
+  G4int pneut = aFragment.GetNumberOfParticles()-pplus;
+  return (pneut >= (GetA()-GetZ()) && pplus >= GetZ());
+}
 
 G4double G4PreCompoundNucleon::
@@ -36 +57 @@
   if ( !IsItPossible(aFragment) ) return 0.0;
 
-  const G4double r0 = G4PreCompoundParameters::GetAddress()->Getr0();
-
   G4double U = aFragment.GetExcitationEnergy();
   G4double P = aFragment.GetNumberOfParticles();
@@ -43 +62 @@
   G4double N = P + H;
   
-  // g = (6.0/pi2)*a*A
-  //  G4EvaporationLevelDensityParameter theLDP;
   G4double g0 = (6.0/pi2)*aFragment.GetA() * 
     G4PreCompoundParameters::GetAddress()->GetLevelDensity();
-  //    theLDP.LevelDensityParameter(static_cast<G4int>(aFragment.GetA()),static_cast<G4int>(aFragment.GetZ()),U);
+ 
   G4double g1 = (6.0/pi2)*GetRestA() * 
     G4PreCompoundParameters::GetAddress()->GetLevelDensity();
-    //    theLDP.LevelDensityParameter(static_cast<G4int>(GetRestA()),static_cast<G4int>(GetRestZ()),U);
 
   G4double A0 = ((P*P+H*H+P-H)/4.0 - H/2.0)/g0;
-  G4double A1 = (A0*g0 - P/2.0)/g1;
+
+  G4double A1 = (A0 - P/2.0)/g1;
   
   G4double E0 = std::max(0.0,U - A0);
@@ -61 +78 @@
 
 
-  G4double Probability = 2.0/(hbarc*hbarc*hbarc) * GetReducedMass() * 
-      GetRj(aFragment.GetNumberOfParticles(), aFragment.GetNumberOfCharged()) * r0 * r0 * std::pow(GetRestA(),2.0/3.0) * GetAlpha() * (eKin + GetBeta()) *
-      P*(N-1.0) * std::pow(g1*E1/(g0*E0),N-2.0)/E0 *
-      g1/(g0*g0);
+  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());
+  }
 
   return Probability;
 }
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+