G4DeuteronEvaporationProbability.cc

Go to the documentation of this file.
00001 //
00002 // ********************************************************************
00003 // * License and Disclaimer                                           *
00004 // *                                                                  *
00005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
00006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
00007 // * conditions of the Geant4 Software License,  included in the file *
00008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
00009 // * include a list of copyright holders.                             *
00010 // *                                                                  *
00011 // * Neither the authors of this software system, nor their employing *
00012 // * institutes,nor the agencies providing financial support for this *
00013 // * work  make  any representation or  warranty, express or implied, *
00014 // * regarding  this  software system or assume any liability for its *
00015 // * use.  Please see the license in the file  LICENSE  and URL above *
00016 // * for the full disclaimer and the limitation of liability.         *
00017 // *                                                                  *
00018 // * This  code  implementation is the result of  the  scientific and *
00019 // * technical work of the GEANT4 collaboration.                      *
00020 // * By using,  copying,  modifying or  distributing the software (or *
00021 // * any work based  on the software)  you  agree  to acknowledge its *
00022 // * use  in  resulting  scientific  publications,  and indicate your *
00023 // * acceptance of all terms of the Geant4 Software license.          *
00024 // ********************************************************************
00025 //
00026 // $Id$
00027 //
00028 // J.M. Quesada (August2008). Based on:
00029 //
00030 // Hadronic Process: Nuclear De-excitations
00031 // by V. Lara (Oct 1998)
00032 //
00033 // Modified:
00034 // 03-09-2008 J.M. Quesada for external choice of inverse cross section option
00035 // 17-11-2010 V.Ivanchenko integer Z and A
00036 
00037 
00038 #include "G4DeuteronEvaporationProbability.hh"
00039 #include "G4SystemOfUnits.hh"
00040 
00041 G4DeuteronEvaporationProbability::G4DeuteronEvaporationProbability() :
00042     G4EvaporationProbability(2,1,3,&theCoulombBarrier) // A,Z,Gamma (fixed JMQ)
00043 {
00044   ResidualA = ResidualZ = theA = theZ = FragmentA = 0;
00045   ResidualAthrd = FragmentAthrd = 0.0;
00046 }
00047 
00048 G4DeuteronEvaporationProbability::~G4DeuteronEvaporationProbability()
00049 {}
00050 
00051 G4double G4DeuteronEvaporationProbability::CalcAlphaParam(const G4Fragment & fragment) 
00052 {
00053   return 1.0 + CCoeficient(fragment.GetZ_asInt()-GetZ());
00054 }
00055 
00056 G4double G4DeuteronEvaporationProbability::CalcBetaParam(const G4Fragment & ) 
00057 {
00058   return 0.0;
00059 }
00060 
00061 G4double G4DeuteronEvaporationProbability::CCoeficient(G4int aZ) 
00062 {
00063   // Data comes from 
00064   // Dostrovsky, Fraenkel and Friedlander
00065   // Physical Review, vol 116, num. 3 1959
00066   // 
00067   // const G4int size = 5;
00068   // G4double Zlist[5] = { 10.0, 20.0, 30.0, 50.0, 70.0};
00069   // G4double Cp[5] = { 0.50, 0.28, 0.20, 0.15, 0.10};
00070   // C for deuteron is equal to C for protons divided by 2
00071   G4double C = 0.0;
00072         
00073   if (aZ >= 70) {
00074     C = 0.10;
00075   } else {
00076     C = ((((0.15417e-06*aZ) - 0.29875e-04)*aZ + 0.21071e-02)*aZ - 0.66612e-01)*aZ + 0.98375;
00077   }
00078         
00079   return C/2.0;
00080         
00081 }
00082 
00083 
00085 //J. M. Quesada (Dec 2007-June 2008): New inverse reaction cross sections 
00086 //OPT=0 Dostrovski's parameterization
00087 //OPT=1,2 Chatterjee's paramaterization 
00088 //OPT=3,4 Kalbach's parameterization 
00089 // 
00090 G4double 
00091 G4DeuteronEvaporationProbability::CrossSection(const  G4Fragment & fragment, G4double K)
00092 {
00093   theA=GetA();
00094   theZ=GetZ();
00095   ResidualA=fragment.GetA_asInt()-theA;
00096   ResidualZ=fragment.GetZ_asInt()-theZ; 
00097   
00098   ResidualAthrd=fG4pow->Z13(ResidualA);
00099   FragmentA=fragment.GetA_asInt();
00100   FragmentAthrd=fG4pow->Z13(FragmentA);
00101 
00102   if (OPTxs==0) {std::ostringstream errOs;
00103   errOs << "We should'n be here (OPT =0) at evaporation cross section calculation (deuterons)!!"  
00104         <<G4endl;
00105   throw G4HadronicException(__FILE__, __LINE__, errOs.str());
00106   return 0.;}
00107   if( OPTxs==1 || OPTxs==2) return G4DeuteronEvaporationProbability::GetOpt12( K);
00108   else if (OPTxs==3 || OPTxs==4)  return G4DeuteronEvaporationProbability::GetOpt34( K);
00109   else{
00110     std::ostringstream errOs;
00111     errOs << "BAD Deuteron CROSS SECTION OPTION AT EVAPORATION!!"  <<G4endl;
00112     throw G4HadronicException(__FILE__, __LINE__, errOs.str());
00113     return 0.;
00114   }
00115 }
00116 
00117 //
00118 //********************* OPT=1,2 : Chatterjee's cross section ********************
00119 //(fitting to cross section from Bechetti & Greenles OM potential)
00120 
00121 G4double G4DeuteronEvaporationProbability::GetOpt12(G4double K)
00122 {
00123   G4double Kc = K;
00124 
00125   // JMQ xsec is set constat above limit of validity
00126   if (K > 50*MeV) { Kc = 50*MeV; }
00127 
00128   G4double landa ,mu ,nu ,p , Ec,q,r,ji,xs;
00129  
00130   G4double    p0 = -38.21;
00131   G4double    p1 = 922.6;
00132   G4double    p2 = -2804.;
00133   G4double    landa0 = -0.0323;
00134   G4double    landa1 = -5.48;
00135   G4double    mum0 = 336.1;
00136   G4double    mu1 = 0.48;
00137   G4double    nu0 = 524.3;
00138   G4double    nu1 = -371.8;
00139   G4double    nu2 = -5.924;  
00140   G4double    delta=1.2;            
00141 
00142   Ec = 1.44*theZ*ResidualZ/(1.5*ResidualAthrd+delta);
00143   p = p0 + p1/Ec + p2/(Ec*Ec);
00144   landa = landa0*ResidualA + landa1;
00145   G4double resmu1 = fG4pow->powZ(ResidualA,mu1); 
00146   mu = mum0*resmu1;
00147   nu = resmu1*(nu0 + nu1*Ec + nu2*(Ec*Ec));
00148   q = landa - nu/(Ec*Ec) - 2*p*Ec;
00149   r = mu + 2*nu/Ec + p*(Ec*Ec);
00150 
00151   ji=std::max(Kc,Ec);
00152   if(Kc < Ec) { xs = p*Kc*Kc + q*Kc + r;}
00153   else {xs = p*(Kc - ji)*(Kc - ji) + landa*Kc + mu + nu*(2 - Kc/ji)/ji ;}
00154                  
00155   if (xs <0.0) {xs=0.0;}
00156               
00157   return xs;
00158 }
00159 
00160 // *********** OPT=3,4 : Kalbach's cross sections (from PRECO code)*************
00161 G4double G4DeuteronEvaporationProbability::GetOpt34(G4double K)
00162 //     ** d from o.m. of perey and perey
00163 {
00164 
00165   G4double landa, mu, nu, p ,signor(1.),sig;
00166   G4double ec,ecsq,xnulam,etest(0.),a; 
00167   G4double b,ecut,cut,ecut2,geom,elab;
00168 
00169   G4double     flow = 1.e-18;
00170   G4double     spill= 1.e+18;
00171 
00172   G4double     p0 = 0.798;
00173   G4double     p1 = 420.3;
00174   G4double     p2 = -1651.;
00175   G4double     landa0 = 0.00619;
00176   G4double     landa1 = -7.54;
00177   G4double     mum0 = 583.5;
00178   G4double     mu1 = 0.337;
00179   G4double     nu0 = 421.8;
00180   G4double     nu1 = -474.5;
00181   G4double     nu2 = -3.592;      
00182  
00183   G4double     ra=0.80;
00184         
00185   //JMQ 13/02/09 increase of reduced radius to lower the barrier
00186   // ec = 1.44 * theZ * ResidualZ / (1.5*ResidualAthrd+ra);
00187   ec = 1.44 * theZ * ResidualZ / (1.7*ResidualAthrd+ra);
00188   ecsq = ec * ec;
00189   p = p0 + p1/ec + p2/ecsq;
00190   landa = landa0*ResidualA + landa1;
00191   a = fG4pow->powZ(ResidualA,mu1);
00192   mu = mum0 * a;
00193   nu = a* (nu0+nu1*ec+nu2*ecsq);  
00194   xnulam = nu / landa;
00195   if (xnulam > spill) { xnulam=0.; }
00196   if (xnulam >= flow) { etest = 1.2 *std::sqrt(xnulam); }
00197 
00198   a = -2.*p*ec + landa - nu/ecsq;
00199   b = p*ecsq + mu + 2.*nu/ec;
00200   ecut = 0.;
00201   cut = a*a - 4.*p*b;
00202   if (cut > 0.) { ecut = std::sqrt(cut); }
00203   ecut = (ecut-a) / (p+p);
00204   ecut2 = ecut;
00205   //JMQ 290310 for avoiding unphysical increase below minimum (at ecut)
00206   //ecut<0 means that there is no cut with energy axis, i.e. xs is set 
00207   //to 0 bellow minimum
00208   //  if (cut < 0.) ecut2 = ecut - 2.;
00209   if (cut < 0.) { ecut2 = ecut; }
00210   elab = K * FragmentA / G4double(ResidualA);
00211   sig = 0.;
00212 
00213   if (elab <= ec) { //start for E<Ec
00214     if (elab > ecut2) { sig = (p*elab*elab+a*elab+b) * signor; }
00215   }           //end for E<Ec
00216   else {           //start for E>Ec
00217     sig = (landa*elab+mu+nu/elab) * signor;
00218     geom = 0.;
00219     if (xnulam < flow || elab < etest) { return sig; }
00220     geom = std::sqrt(theA*K);
00221     geom = 1.23*ResidualAthrd + ra + 4.573/geom;
00222     geom = 31.416 * geom * geom;
00223     sig = std::max(geom,sig);
00224   }           //end for E>Ec
00225   return sig;
00226 }
00227 

Generated on Mon May 27 17:48:00 2013 for Geant4 by  doxygen 1.4.7