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

Timestamp:

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

Author:

garnier

Message:

update processes

File:

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

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

@@ -24 +24 @@
 // ********************************************************************
 //
-//
-// $Id: G4AlphaEvaporationProbability.cc,v 1.4 2006/06/29 20:10:19 gunter Exp $
-// GEANT4 tag $Name: geant4-09-01-patch-02 $
+//J.M. Quesada (August2008). Based on:
 //
 // Hadronic Process: Nuclear De-excitations
-// by V. Lara (Nov 1999)
-//
-
+// by V. Lara (Oct 1998)
+//
+// Modif (03 September 2008) by J. M. Quesada for external choice of inverse 
+// cross section option
 
 #include "G4AlphaEvaporationProbability.hh"
 
 G4AlphaEvaporationProbability::G4AlphaEvaporationProbability() :
-    G4EvaporationProbability(4,2,4) // A,Z,Gamma
-{
-    //  const G4int NumExcitedStates = 31+1;
-    std::vector<G4double>::size_type NumExcitedStatesEnergy = 31+1;
-    std::vector<G4int>::size_type NumExcitedStatesSpin = 31+1;
-    ExcitEnergies.reserve(NumExcitedStatesEnergy);
-    ExcitSpins.reserve(NumExcitedStatesSpin);
-    ExcitEnergies.insert(ExcitEnergies.begin(),NumExcitedStatesEnergy,0.0);
-    ExcitSpins.insert(ExcitSpins.begin(),NumExcitedStatesSpin,0);
-//    for (G4int i = 0; i < NumExcitedStates; i++) {
-//      ExcitEnergies(i) = 0.0;
-//      ExcitSpins(i) = 0;
-//    }
-
-    ExcitEnergies[18] = 7.98*MeV;
-    ExcitEnergies[20] = 6.90*MeV;
-    ExcitEnergies[25] = 5.83*MeV;
-    ExcitEnergies[26] = 8.57*MeV;
-    ExcitEnergies[31] = 5.33*MeV;
-
-    ExcitSpins[18] = 4;
-    ExcitSpins[20] = 6;
-    ExcitSpins[25] = 7;
-    ExcitSpins[26] = 4;
-    ExcitSpins[31] = 13;
-        
-    SetExcitationEnergiesPtr(&ExcitEnergies);
-    SetExcitationSpinsPtr(&ExcitSpins);         
+    G4EvaporationProbability(4,2,1,&theCoulombBarrier) // A,Z,Gamma,&theCoumlombBarrier
+{
+        
 }
 
@@ -93 +67 @@
 }
 
-G4double G4AlphaEvaporationProbability::CCoeficient(const G4double aZ) const
+
+ 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(const G4double aZ) 
 {
     // Data comes from 
@@ -116 +97 @@
     return C;
 }
+
+///////////////////////////////////////////////////////////////////////////////////
+//J. M. Quesada (Dec 2007-June 2008): New inverse reaction cross sections 
+//OPT=0 Dostrovski's parameterization
+//OPT=1,2 Chatterjee's paramaterization 
+//OPT=3,4 Kalbach's parameterization 
+// 
+G4double G4AlphaEvaporationProbability::CrossSection(const  G4Fragment & fragment, 
+                                                     const  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);
+  
+  
+  if (OPTxs==0) {std::ostringstream errOs;
+    errOs << "We should'n be here (OPT =0) at evaporation cross section calculation (Alpha's)!!"  
+          <<G4endl;
+    throw G4HadronicException(__FILE__, __LINE__, errOs.str());
+    return 0.;}
+  
+  if( OPTxs==1 || OPTxs==2) return G4AlphaEvaporationProbability::GetOpt12( K);
+  else if (OPTxs==3 || OPTxs==4)  return G4AlphaEvaporationProbability::GetOpt34( K);
+  else{
+    std::ostringstream errOs;
+    errOs << "BAD Alpha CROSS SECTION OPTION AT EVAPORATION!!"  <<G4endl;
+    throw G4HadronicException(__FILE__, __LINE__, errOs.str());
+    return 0.;
+  }
+}
+
+
+//
+//********************* 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 Kc=K;
+  
+  // JMQ xsec is set constat above limit of validity
+  if (K>50) Kc=50;
+  
+  G4double landa ,mu ,nu ,p , Ec,q,r,ji,xs;
+  
+  G4double     p0 = 10.95;
+  G4double     p1 = -85.2;
+  G4double     p2 = 1146.;
+  G4double     landa0 = 0.0643;
+  G4double     landa1 = -13.96;
+  G4double     mu0 = 781.2;
+  G4double     mu1 = 0.29;
+  G4double     nu0 = -304.7;
+  G4double     nu1 = -470.0;
+  G4double     nu2 = -8.580;   
+  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));
+  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)
+// 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     p1 = -85.2;
+  G4double     p2 = 1146.;
+  G4double     landa0 = 0.0643;
+  G4double     landa1 = -13.96;
+  G4double     mu0 = 781.2;
+  G4double     mu1 = 0.29;
+  G4double     nu0 = -304.7;
+  G4double     nu1 = -470.0;
+  G4double     nu2 = -8.580;        
+  
+  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);
+  ec = 1.44 * theZ * ResidualZ / (1.7*ResidualAthrd+ra);
+  ecsq = ec * ec;
+  p = p0 + p1/ec + p2/ecsq;
+  landa = landa0*ResidualA + landa1;
+  a = std::pow(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);
+  
+  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);
+  ecut = (ecut-a) / (p+p);
+  ecut2 = ecut;
+  if (cut < 0.) ecut2 = ecut - 2.;
+  elab = K * FragmentA / ResidualA;
+  sig = 0.;
+  
+  if (elab <= ec) { //start for E<Ec
+    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;
+    geom = std::sqrt(theA*K);
+    geom = 1.23*ResidualAthrd + ra + 4.573/geom;
+    geom = 31.416 * geom * geom;
+    sig = std::max(geom,sig);
+  }           //end for E>Ec
+  return sig;
+  
+}
+
+//   ************************** end of cross sections *******************************