Changeset 1347 for trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4AlphaEvaporationProbability.cc

Timestamp:

Dec 22, 2010, 3:52:27 PM (13 years ago)

Author:

garnier

Message:

geant4 tag 9.4

File:

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

trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4AlphaEvaporationProbability.cc

@@ -24 +24 @@
 // ********************************************************************
 //
-//J.M. Quesada (August2008). Based on:
+// $Id: G4AlphaEvaporationProbability.cc,v 1.18 2010/11/17 11:06:03 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-04-ref-00 $
+//
+// J.M. Quesada (August2008). Based on:
 //
 // Hadronic Process: Nuclear De-excitations
 // by V. Lara (Oct 1998)
 //
-// Modif (03 September 2008) by J. M. Quesada for external choice of inverse 
-// cross section option
+// Modified:
+// 03-09-2008 J.M. Quesada for external choice of inverse cross section option
+// 17-11-2010 V.Ivanchenko integer Z and A
 
 #include "G4AlphaEvaporationProbability.hh"
@@ -36 +40 @@
 G4AlphaEvaporationProbability::G4AlphaEvaporationProbability() :
     G4EvaporationProbability(4,2,1,&theCoulombBarrier) // A,Z,Gamma,&theCoumlombBarrier
-{
+{}
+
+G4AlphaEvaporationProbability::~G4AlphaEvaporationProbability()
+{}
+
+G4double G4AlphaEvaporationProbability::CalcAlphaParam(const G4Fragment & fragment)  
+  { return 1.0 + CCoeficient(fragment.GetZ_asInt()-GetZ());}
         
-}
-
-
-G4AlphaEvaporationProbability::G4AlphaEvaporationProbability(const G4AlphaEvaporationProbability &): G4EvaporationProbability()
-{
-    throw G4HadronicException(__FILE__, __LINE__, "G4AlphaEvaporationProbability::copy_constructor meant to not be accessable");
-}
-
-
-
-
-const G4AlphaEvaporationProbability & G4AlphaEvaporationProbability::
-operator=(const G4AlphaEvaporationProbability &) 
-{
-    throw G4HadronicException(__FILE__, __LINE__, "G4AlphaEvaporationProbability::operator= meant to not be accessable");
-    return *this;
-}
-
-
-G4bool G4AlphaEvaporationProbability::operator==(const G4AlphaEvaporationProbability &) const
-{
-    return false;
-}
-
-G4bool G4AlphaEvaporationProbability::operator!=(const G4AlphaEvaporationProbability &) const
-{
-    return true;
-}
-
-
- G4double G4AlphaEvaporationProbability::CalcAlphaParam(const G4Fragment & fragment)  
-  { return 1.0 + CCoeficient(static_cast<G4double>(fragment.GetZ()-GetZ()));}
+G4double G4AlphaEvaporationProbability::CalcBetaParam(const G4Fragment &) 
+  { return 0.0; }
+
+G4double G4AlphaEvaporationProbability::CCoeficient(G4int aZ) 
+{
+  // Data comes from 
+  // Dostrovsky, Fraenkel and Friedlander
+  // Physical Review, vol 116, num. 3 1959
+  // 
+  // const G4int size = 5;
+  // G4double Zlist[5] = { 10.0, 20.0, 30.0, 50.0, 70.0};
+  //    G4double Calpha[5] = { 0.10, 0.10, 0.10, 0.08, 0.06};
+  G4double C = 0.0;
         
- G4double G4AlphaEvaporationProbability::CalcBetaParam(const G4Fragment &) 
-  { return 0.0; }
-
-G4double G4AlphaEvaporationProbability::CCoeficient(const G4double aZ) 
-{
-    // Data comes from 
-    // Dostrovsky, Fraenkel and Friedlander
-    // Physical Review, vol 116, num. 3 1959
-    // 
-    // const G4int size = 5;
-    // G4double Zlist[5] = { 10.0, 20.0, 30.0, 50.0, 70.0};
-    //  G4double Calpha[5] = { 0.10, 0.10, 0.10, 0.08, 0.06};
-    G4double C = 0.0;
-        
-        
-    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 C;
+  if (aZ <= 30) 
+    {
+      C = 0.10;
+    }
+  else if (aZ <= 50)
+    {
+      C = 0.1 - (aZ-30)*0.001;
+    }
+  else if (aZ < 70)
+    {
+      C = 0.08 - (aZ-50)*0.001;
+    }
+  else 
+    {
+      C = 0.06;
+    }
+  return C;
 }
 
@@ -104 +87 @@
 //OPT=3,4 Kalbach's parameterization 
 // 
-G4double G4AlphaEvaporationProbability::CrossSection(const  G4Fragment & fragment, 
-                                                     const  G4double K)
+G4double 
+G4AlphaEvaporationProbability::CrossSection(const G4Fragment & fragment, G4double K)
 {
   theA=GetA();
   theZ=GetZ();
-  ResidualA=fragment.GetA()-theA;
-  ResidualZ=fragment.GetZ()-theZ; 
-  
-  ResidualAthrd=std::pow(ResidualA,0.33333);
-  FragmentA=fragment.GetA();
-  FragmentAthrd=std::pow(FragmentA,0.33333);
-  
+  ResidualA=fragment.GetA_asInt()-theA;
+  ResidualZ=fragment.GetZ_asInt()-theZ; 
+  
+  ResidualAthrd=fG4pow->Z13(ResidualA);
+  FragmentA=fragment.GetA_asInt();
+  FragmentAthrd=fG4pow->Z13(FragmentA);
   
   if (OPTxs==0) {std::ostringstream errOs;
@@ -133 +115 @@
 }
 
-
-//
-//********************* OPT=1,2 : Chatterjee's cross section ************************ 
+//
+//********************* OPT=1,2 : Chatterjee's cross section ********************
 //(fitting to cross section from Bechetti & Greenles OM potential)
 
-G4double G4AlphaEvaporationProbability::GetOpt12(const  G4double K)
-{
-// c     ** alpha from huizenga and igo
+G4double G4AlphaEvaporationProbability::GetOpt12(G4double K)
+{
   G4double Kc=K;
-  
-  // JMQ xsec is set constat above limit of validity
-  if (K>50) Kc=50;
-  
+
+  // JMQ xsec is set constant above limit of validity
+  if (K > 50*MeV) { Kc = 50*MeV; }
+
   G4double landa ,mu ,nu ,p , Ec,q,r,ji,xs;
-  
+
   G4double     p0 = 10.95;
   G4double     p1 = -85.2;
@@ -160 +140 @@
   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 = fG4pow->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);
-  
+
   ji=std::max(Kc,Ec);
   if(Kc < Ec) { xs = p*Kc*Kc + q*Kc + r;}
   else {xs = p*(Kc - ji)*(Kc - ji) + landa*Kc + mu + nu*(2 - Kc/ji)/ji ;}
-                 
+  
   if (xs <0.0) {xs=0.0;}
-  
+              
   return xs;
-  
 }
 
 // *********** OPT=3,4 : Kalbach's cross sections (from PRECO code)*************
-G4double G4AlphaEvaporationProbability::GetOpt34(const  G4double K)
+G4double G4AlphaEvaporationProbability::GetOpt34(G4double K)
 // c     ** alpha from huizenga and igo
 {
-  
   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 = 10.95;
+  G4double     p0 = 10.95;
   G4double     p1 = -85.2;
   G4double     p2 = 1146.;
@@ -204 +181 @@
   
   G4double      ra=1.20;
-  
+        
   //JMQ 13/02/09 increase of reduced radius to lower the barrier
   // ec = 1.44 * theZ * ResidualZ / (1.5*ResidualAthrd+ra);
@@ -211 +188 @@
   p = p0 + p1/ec + p2/ecsq;
   landa = landa0*ResidualA + landa1;
-  a = std::pow(ResidualA,mu1);
+  a = fG4pow->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;
   b = p*ecsq + mu + 2.*nu/ec;
   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;
@@ -245 +223 @@
   }           //end for E>Ec
   return sig;
-  
-}
-
-//   ************************** end of cross sections ******************************* 
+}
+