Ignore:
Timestamp:
Dec 22, 2010, 3:52:27 PM (13 years ago)
Author:
garnier
Message:

geant4 tag 9.4

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4He3EvaporationProbability.cc

    r1315 r1347  
    2424// ********************************************************************
    2525//
    26 //J.M. Quesada (August2008). Based on:
     26// $Id: G4He3EvaporationProbability.cc,v 1.18 2010/11/17 11:06:03 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-04-ref-00 $
     28//
     29// J.M. Quesada (August2008). Based on:
    2730//
    2831// Hadronic Process: Nuclear De-excitations
    2932// by V. Lara (Oct 1998)
    3033//
    31 // Modif (03 September 2008) by J. M. Quesada for external choice of inverse
    32 // cross section option
     34// Modified:
     35// 03-09-2008 J.M. Quesada for external choice of inverse cross section option
     36// 17-11-2010 V.Ivanchenko integer Z and A
    3337
    3438#include "G4He3EvaporationProbability.hh"
    … …  
    3640G4He3EvaporationProbability::G4He3EvaporationProbability() :
    3741   G4EvaporationProbability(3,2,2,&theCoulombBarrier) // A,Z,Gamma,&theCoulombBarrier
    38 {
    39 
    40 }
    41 
    42 G4He3EvaporationProbability::G4He3EvaporationProbability(const G4He3EvaporationProbability &) : G4EvaporationProbability()
    43 {
    44     throw G4HadronicException(__FILE__, __LINE__, "G4He3EvaporationProbability::copy_constructor meant to not be accessable");
    45 }
    46 
    47 
    48 
    49 
    50 const G4He3EvaporationProbability & G4He3EvaporationProbability::
    51 operator=(const G4He3EvaporationProbability &)
    52 {
    53     throw G4HadronicException(__FILE__, __LINE__, "G4He3EvaporationProbability::operator= meant to not be accessable");
    54     return *this;
    55 }
    56 
    57 
    58 G4bool G4He3EvaporationProbability::operator==(const G4He3EvaporationProbability &) const
    59 {
    60     return false;
    61 }
    62 
    63 G4bool G4He3EvaporationProbability::operator!=(const G4He3EvaporationProbability &) const
    64 {
    65     return true;
    66 }
    67 
    68   G4double G4He3EvaporationProbability::CalcAlphaParam(const G4Fragment & fragment) 
    69   { return 1.0 + CCoeficient(static_cast<G4double>(fragment.GetZ()-GetZ()));}
     42{}
     43
     44G4He3EvaporationProbability::~G4He3EvaporationProbability()
     45{}
     46
     47G4double G4He3EvaporationProbability::CalcAlphaParam(const G4Fragment & fragment) 
     48  { return 1.0 + CCoeficient(fragment.GetZ_asInt()-GetZ());}
    7049       
    71   G4double G4He3EvaporationProbability::CalcBetaParam(const G4Fragment & ) 
     50G4double G4He3EvaporationProbability::CalcBetaParam(const G4Fragment & ) 
    7251  { return 0.0; }
    7352
    7453
    75 G4double G4He3EvaporationProbability::CCoeficient(const G4double aZ)
    76 {
    77     // Data comes from
    78     // Dostrovsky, Fraenkel and Friedlander
    79     // Physical Review, vol 116, num. 3 1959
    80     //
    81     // const G4int size = 5;
    82     // G4double Zlist[5] = { 10.0, 20.0, 30.0, 50.0, 70.0};
    83     //  G4double Calpha[5] = { 0.10, 0.10, 0.10, 0.08, 0.06};
    84     // C for He3 is equal to C for alpha times 4/3
    85     G4double C = 0.0;
     54G4double G4He3EvaporationProbability::CCoeficient(G4int aZ)
     55{
     56  // Data comes from
     57  // Dostrovsky, Fraenkel and Friedlander
     58  // Physical Review, vol 116, num. 3 1959
     59  //
     60  // const G4int size = 5;
     61  // G4double Zlist[5] = { 10.0, 20.0, 30.0, 50.0, 70.0};
     62  //    G4double Calpha[5] = { 0.10, 0.10, 0.10, 0.08, 0.06};
     63  // C for He3 is equal to C for alpha times 4/3
     64  G4double C = 0.0;
    8665       
    87        
    88     if (aZ <= 30) {
    89         C = 0.10;
    90     } else if (aZ <= 50) {
    91         C = 0.1 + -((aZ-50.)/20.)*0.02;
    92     } else if (aZ < 70) {
    93         C = 0.08 + -((aZ-70.)/20.)*0.02;
    94     } else {
    95         C = 0.06;
    96     }
    97     return C*(4.0/3.0);
     66  if (aZ <= 30)
     67    {
     68      C = 0.10;
     69    }
     70  else if (aZ <= 50)
     71    {
     72      C = 0.1 - (aZ - 30)*0.001;
     73    }
     74  else if (aZ < 70)
     75    {
     76      C = 0.08 - (aZ - 50)*0.001;
     77    }
     78  else
     79    {
     80      C = 0.06;
     81    }
     82  return C*(4.0/3.0);
    9883}
    9984
    … …  
    10489//OPT=3,4 Kalbach's parameterization
    10590//
    106 G4double G4He3EvaporationProbability::CrossSection(const  G4Fragment & fragment, const  G4double K)
    107 {
    108        theA=GetA();
    109        theZ=GetZ();
    110        ResidualA=fragment.GetA()-theA;
    111        ResidualZ=fragment.GetZ()-theZ;
    112  
    113        ResidualAthrd=std::pow(ResidualA,0.33333);
    114        FragmentA=fragment.GetA();
    115        FragmentAthrd=std::pow(FragmentA,0.33333);
    116 
    117 
    118        if (OPTxs==0) {std::ostringstream errOs;
    119          errOs << "We should'n be here (OPT =0) at evaporation cross section calculation (He3's)!!" 
    120                <<G4endl;
    121          throw G4HadronicException(__FILE__, __LINE__, errOs.str());
    122          return 0.;}
    123        if( OPTxs==1 || OPTxs==2) return G4He3EvaporationProbability::GetOpt12( K);
    124        else if (OPTxs==3 || OPTxs==4)  return G4He3EvaporationProbability::GetOpt34( K);
    125        else{
    126          std::ostringstream errOs;
    127          errOs << "BAD He3's CROSS SECTION OPTION AT EVAPORATION!!"  <<G4endl;
    128          throw G4HadronicException(__FILE__, __LINE__, errOs.str());
    129          return 0.;
    130        }
    131 }
    132 //
    133 //********************* OPT=1,2 : Chatterjee's cross section ************************
     91G4double
     92G4He3EvaporationProbability::CrossSection(const  G4Fragment & fragment, G4double K)
     93{
     94
     95  theA=GetA();
     96  theZ=GetZ();
     97  ResidualA=fragment.GetA_asInt()-theA;
     98  ResidualZ=fragment.GetZ_asInt()-theZ;
     99 
     100  ResidualAthrd=fG4pow->Z13(ResidualA);
     101  FragmentA=fragment.GetA_asInt();
     102  FragmentAthrd=fG4pow->Z13(FragmentA);
     103
     104  if (OPTxs==0) {std::ostringstream errOs;
     105  errOs << "We should'n be here (OPT =0) at evaporation cross section calculation (He3's)!!" 
     106        <<G4endl;
     107  throw G4HadronicException(__FILE__, __LINE__, errOs.str());
     108  return 0.;}
     109  if( OPTxs==1 || OPTxs==2) return G4He3EvaporationProbability::GetOpt12( K);
     110  else if (OPTxs==3 || OPTxs==4)  return G4He3EvaporationProbability::GetOpt34( K);
     111  else{
     112    std::ostringstream errOs;
     113    errOs << "BAD He3's CROSS SECTION OPTION AT EVAPORATION!!"  <<G4endl;
     114    throw G4HadronicException(__FILE__, __LINE__, errOs.str());
     115    return 0.;
     116  }
     117}
     118
     119//********************* OPT=1,2 : Chatterjee's cross section *****************
    134120//(fitting to cross section from Bechetti & Greenles OM potential)
    135121
    136122G4double G4He3EvaporationProbability::GetOpt12(const  G4double K)
    137123{
    138 
    139 G4double Kc=K;
    140 
    141 // JMQ xsec is set constat above limit of validity
    142  if (K>50) Kc=50;
    143 
    144  G4double landa ,mu ,nu ,p , Ec,q,r,ji,xs;
    145 
    146  G4double     p0 = -3.06;
    147  G4double     p1 = 278.5;
    148  G4double     p2 = -1389.;
    149  G4double     landa0 = -0.00535;
    150  G4double     landa1 = -11.16;
    151  G4double     mu0 = 555.5;
    152  G4double     mu1 = 0.40;
    153  G4double     nu0 = 687.4;
    154  G4double     nu1 = -476.3;
    155  G4double     nu2 = 0.509;   
    156  G4double     delta=1.2;             
    157 
    158       Ec = 1.44*theZ*ResidualZ/(1.5*ResidualAthrd+delta);
    159       p = p0 + p1/Ec + p2/(Ec*Ec);
    160       landa = landa0*ResidualA + landa1;
    161       mu = mu0*std::pow(ResidualA,mu1);
    162       nu = std::pow(ResidualA,mu1)*(nu0 + nu1*Ec + nu2*(Ec*Ec));
    163       q = landa - nu/(Ec*Ec) - 2*p*Ec;
    164       r = mu + 2*nu/Ec + p*(Ec*Ec);
    165 
    166    ji=std::max(Kc,Ec);
    167    if(Kc < Ec) { xs = p*Kc*Kc + q*Kc + r;}
    168    else {xs = p*(Kc - ji)*(Kc - ji) + landa*Kc + mu + nu*(2 - Kc/ji)/ji ;}
    169    
    170    if (xs <0.0) {xs=0.0;}
     124  G4double Kc = K;
     125
     126  // JMQ xsec is set constat above limit of validity
     127  if (K > 50*MeV) { Kc = 50*MeV; }
     128
     129  G4double landa ,mu ,nu ,p , Ec,q,r,ji,xs;
     130
     131  G4double     p0 = -3.06;
     132  G4double     p1 = 278.5;
     133  G4double     p2 = -1389.;
     134  G4double     landa0 = -0.00535;
     135  G4double     landa1 = -11.16;
     136  G4double     mu0 = 555.5;
     137  G4double     mu1 = 0.40;
     138  G4double     nu0 = 687.4;
     139  G4double     nu1 = -476.3;
     140  G4double     nu2 = 0.509;   
     141  G4double     delta=1.2;             
     142
     143  Ec = 1.44*theZ*ResidualZ/(1.5*ResidualAthrd+delta);
     144  p = p0 + p1/Ec + p2/(Ec*Ec);
     145  landa = landa0*ResidualA + landa1;
     146
     147  G4double resmu1 = fG4pow->powZ(ResidualA,mu1);
     148  mu = mu0*resmu1;
     149  nu = resmu1*(nu0 + nu1*Ec + nu2*(Ec*Ec));
     150  q = landa - nu/(Ec*Ec) - 2*p*Ec;
     151  r = mu + 2*nu/Ec + p*(Ec*Ec);
     152 
     153  ji=std::max(Kc,Ec);
     154  if(Kc < Ec) { xs = p*Kc*Kc + q*Kc + r;}
     155  else {xs = p*(Kc - ji)*(Kc - ji) + landa*Kc + mu + nu*(2 - Kc/ji)/ji ;}
     156 
     157  if (xs <0.0) {xs=0.0;}
    171158             
    172    return xs;
     159  return xs;
    173160
    174161}
    … …  
    178165//c     ** 3he from o.m. of gibson et al
    179166{
    180  
    181167  G4double landa, mu, nu, p , signor(1.),sig;
    182168  G4double ec,ecsq,xnulam,etest(0.),a;
    183169  G4double b,ecut,cut,ecut2,geom,elab;
    184  
    185  
     170
    186171  G4double     flow = 1.e-18;
    187172  G4double     spill= 1.e+18;
    188 
    189173
    190174  G4double     p0 = -2.88;
    … …  
    200184 
    201185  G4double      ra=0.80;
    202  
     186       
    203187  //JMQ 13/02/09 increase of reduced radius to lower the barrier
    204188  // ec = 1.44 * theZ * ResidualZ / (1.5*ResidualAthrd+ra);
    … …  
    207191  p = p0 + p1/ec + p2/ecsq;
    208192  landa = landa0*ResidualA + landa1;
    209   a = std::pow(ResidualA,mu1);
     193  a = fG4pow->powZ(ResidualA,mu1);
    210194  mu = mu0 * a;
    211195  nu = a* (nu0+nu1*ec+nu2*ecsq); 
    212196  xnulam = nu / landa;
    213   if (xnulam > spill) xnulam=0.;
    214   if (xnulam >= flow) etest = 1.2 *std::sqrt(xnulam);
     197  if (xnulam > spill) { xnulam=0.; }
     198  if (xnulam >= flow) { etest = 1.2 *std::sqrt(xnulam); }
    215199 
    216200  a = -2.*p*ec + landa - nu/ecsq;
    … …  
    221205  ecut = (ecut-a) / (p+p);
    222206  ecut2 = ecut;
    223 //JMQ 290310 for avoiding unphysical increase below minimum (at ecut)
    224 //ecut<0 means that there is no cut with energy axis, i.e. xs is set to 0 bellow minimum
    225 //  if (cut < 0.) ecut2 = ecut - 2.;
    226   if (cut < 0.) ecut2 = ecut;
    227   elab = K * FragmentA / ResidualA;
     207  //JMQ 290310 for avoiding unphysical increase below minimum (at ecut)
     208  // ecut<0 means that there is no cut with energy axis, i.e. xs is set
     209  // to 0 bellow minimum
     210  //  if (cut < 0.) ecut2 = ecut - 2.;
     211  if (cut < 0.) { ecut2 = ecut; }
     212  elab = K * FragmentA /G4double(ResidualA);
    228213  sig = 0.;
    229 
     214 
    230215  if (elab <= ec) { //start for E<Ec
    231     if (elab > ecut2)  sig = (p*elab*elab+a*elab+b) * signor;
     216    if (elab > ecut2) { sig = (p*elab*elab+a*elab+b) * signor; }
    232217  }           //end for E<Ec
    233218  else {           //start for E>Ec
    234219    sig = (landa*elab+mu+nu/elab) * signor;
    235220    geom = 0.;
    236     if (xnulam < flow || elab < etest) return sig;
     221    if (xnulam < flow || elab < etest) { return sig; }
    237222    geom = std::sqrt(theA*K);
    238223    geom = 1.23*ResidualAthrd + ra + 4.573/geom;
    … …  
    243228 
    244229}
    245 
    246 //   ************************** end of cross sections *******************************
    247 
Note: See TracChangeset for help on using the changeset viewer.