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

Timestamp:

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

Author:

garnier

Message:

update ti head

File:

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

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;
+}