Changeset 1340 for trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/src/G4PreCompoundAlpha.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/G4PreCompoundAlpha.cc (modified) (12 diffs)

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