Ignore:
Timestamp:
Apr 6, 2009, 12:30:29 PM (15 years ago)
Author:
garnier
Message:

update processes

Location:
trunk/source/processes/hadronic/models/de_excitation/evaporation/src
Files:
17 edited

Legend:

Unmodified
Added
Removed

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

    r819 r962  
    2626//
    2727// $Id: G4AlphaEvaporationChannel.cc,v 1.4 2006/06/29 20:10:17 gunter Exp $
    28 // GEANT4 tag $Name: $
     28// GEANT4 tag $Name: geant4-09-02-ref-02 $
    2929//
    3030// Hadronic Process: Nuclear De-excitations

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

    r819 r962  
    2424// ********************************************************************
    2525//
    26 //
    27 // $Id: G4AlphaEvaporationProbability.cc,v 1.4 2006/06/29 20:10:19 gunter Exp $
    28 // GEANT4 tag $Name: geant4-09-01-patch-02 $
     26//J.M. Quesada (August2008). Based on:
    2927//
    3028// Hadronic Process: Nuclear De-excitations
    31 // by V. Lara (Nov 1999)
    32 //
    33 
     29// by V. Lara (Oct 1998)
     30//
     31// Modif (03 September 2008) by J. M. Quesada for external choice of inverse
     32// cross section option
    3433
    3534#include "G4AlphaEvaporationProbability.hh"
    3635
    3736G4AlphaEvaporationProbability::G4AlphaEvaporationProbability() :
    38     G4EvaporationProbability(4,2,4) // A,Z,Gamma
    39 {
    40     //  const G4int NumExcitedStates = 31+1;
    41     std::vector<G4double>::size_type NumExcitedStatesEnergy = 31+1;
    42     std::vector<G4int>::size_type NumExcitedStatesSpin = 31+1;
    43     ExcitEnergies.reserve(NumExcitedStatesEnergy);
    44     ExcitSpins.reserve(NumExcitedStatesSpin);
    45     ExcitEnergies.insert(ExcitEnergies.begin(),NumExcitedStatesEnergy,0.0);
    46     ExcitSpins.insert(ExcitSpins.begin(),NumExcitedStatesSpin,0);
    47 //    for (G4int i = 0; i < NumExcitedStates; i++) {
    48 //      ExcitEnergies(i) = 0.0;
    49 //      ExcitSpins(i) = 0;
    50 //    }
    51 
    52     ExcitEnergies[18] = 7.98*MeV;
    53     ExcitEnergies[20] = 6.90*MeV;
    54     ExcitEnergies[25] = 5.83*MeV;
    55     ExcitEnergies[26] = 8.57*MeV;
    56     ExcitEnergies[31] = 5.33*MeV;
    57 
    58     ExcitSpins[18] = 4;
    59     ExcitSpins[20] = 6;
    60     ExcitSpins[25] = 7;
    61     ExcitSpins[26] = 4;
    62     ExcitSpins[31] = 13;
    63        
    64     SetExcitationEnergiesPtr(&ExcitEnergies);
    65     SetExcitationSpinsPtr(&ExcitSpins);         
     37    G4EvaporationProbability(4,2,1,&theCoulombBarrier) // A,Z,Gamma,&theCoumlombBarrier
     38{
     39       
    6640}
    6741
    … …  
    9367}
    9468
    95 G4double G4AlphaEvaporationProbability::CCoeficient(const G4double aZ) const
     69
     70 G4double G4AlphaEvaporationProbability::CalcAlphaParam(const G4Fragment & fragment) 
     71  { return 1.0 + CCoeficient(static_cast<G4double>(fragment.GetZ()-GetZ()));}
     72       
     73 G4double G4AlphaEvaporationProbability::CalcBetaParam(const G4Fragment &)
     74  { return 0.0; }
     75
     76G4double G4AlphaEvaporationProbability::CCoeficient(const G4double aZ)
    9677{
    9778    // Data comes from
    … …  
    11697    return C;
    11798}
     99
     100///////////////////////////////////////////////////////////////////////////////////
     101//J. M. Quesada (Dec 2007-June 2008): New inverse reaction cross sections
     102//OPT=0 Dostrovski's parameterization
     103//OPT=1,2 Chatterjee's paramaterization
     104//OPT=3,4 Kalbach's parameterization
     105//
     106G4double G4AlphaEvaporationProbability::CrossSection(const  G4Fragment & fragment,
     107                                                     const  G4double K)
     108{
     109  theA=GetA();
     110  theZ=GetZ();
     111  ResidualA=fragment.GetA()-theA;
     112  ResidualZ=fragment.GetZ()-theZ;
     113 
     114  ResidualAthrd=std::pow(ResidualA,0.33333);
     115  FragmentA=fragment.GetA();
     116  FragmentAthrd=std::pow(FragmentA,0.33333);
     117 
     118 
     119  if (OPTxs==0) {std::ostringstream errOs;
     120    errOs << "We should'n be here (OPT =0) at evaporation cross section calculation (Alpha's)!!" 
     121          <<G4endl;
     122    throw G4HadronicException(__FILE__, __LINE__, errOs.str());
     123    return 0.;}
     124 
     125  if( OPTxs==1 || OPTxs==2) return G4AlphaEvaporationProbability::GetOpt12( K);
     126  else if (OPTxs==3 || OPTxs==4)  return G4AlphaEvaporationProbability::GetOpt34( K);
     127  else{
     128    std::ostringstream errOs;
     129    errOs << "BAD Alpha CROSS SECTION OPTION AT EVAPORATION!!"  <<G4endl;
     130    throw G4HadronicException(__FILE__, __LINE__, errOs.str());
     131    return 0.;
     132  }
     133}
     134
     135
     136//
     137//********************* OPT=1,2 : Chatterjee's cross section ************************
     138//(fitting to cross section from Bechetti & Greenles OM potential)
     139
     140G4double G4AlphaEvaporationProbability::GetOpt12(const  G4double K)
     141{
     142// c     ** alpha from huizenga and igo
     143  G4double Kc=K;
     144 
     145  // JMQ xsec is set constat above limit of validity
     146  if (K>50) Kc=50;
     147 
     148  G4double landa ,mu ,nu ,p , Ec,q,r,ji,xs;
     149 
     150  G4double     p0 = 10.95;
     151  G4double     p1 = -85.2;
     152  G4double     p2 = 1146.;
     153  G4double     landa0 = 0.0643;
     154  G4double     landa1 = -13.96;
     155  G4double     mu0 = 781.2;
     156  G4double     mu1 = 0.29;
     157  G4double     nu0 = -304.7;
     158  G4double     nu1 = -470.0;
     159  G4double     nu2 = -8.580;   
     160  G4double     delta=1.2;         
     161
     162 
     163  Ec = 1.44*theZ*ResidualZ/(1.5*ResidualAthrd+delta);
     164  p = p0 + p1/Ec + p2/(Ec*Ec);
     165  landa = landa0*ResidualA + landa1;
     166  mu = mu0*std::pow(ResidualA,mu1);
     167  nu = std::pow(ResidualA,mu1)*(nu0 + nu1*Ec + nu2*(Ec*Ec));
     168  q = landa - nu/(Ec*Ec) - 2*p*Ec;
     169  r = mu + 2*nu/Ec + p*(Ec*Ec);
     170 
     171  ji=std::max(Kc,Ec);
     172  if(Kc < Ec) { xs = p*Kc*Kc + q*Kc + r;}
     173  else {xs = p*(Kc - ji)*(Kc - ji) + landa*Kc + mu + nu*(2 - Kc/ji)/ji ;}
     174                 
     175  if (xs <0.0) {xs=0.0;}
     176 
     177  return xs;
     178 
     179}
     180
     181// *********** OPT=3,4 : Kalbach's cross sections (from PRECO code)*************
     182G4double G4AlphaEvaporationProbability::GetOpt34(const  G4double K)
     183// c     ** alpha from huizenga and igo
     184{
     185 
     186  G4double landa, mu, nu, p , signor(1.),sig;
     187  G4double ec,ecsq,xnulam,etest(0.),a;
     188  G4double b,ecut,cut,ecut2,geom,elab;
     189
     190
     191  G4double     flow = 1.e-18;
     192  G4double     spill= 1.e+18;
     193
     194  G4double       p0 = 10.95;
     195  G4double     p1 = -85.2;
     196  G4double     p2 = 1146.;
     197  G4double     landa0 = 0.0643;
     198  G4double     landa1 = -13.96;
     199  G4double     mu0 = 781.2;
     200  G4double     mu1 = 0.29;
     201  G4double     nu0 = -304.7;
     202  G4double     nu1 = -470.0;
     203  G4double     nu2 = -8.580;       
     204 
     205  G4double      ra=1.20;
     206 
     207  //JMQ 13/02/09 increase of reduced radius to lower the barrier
     208  // ec = 1.44 * theZ * ResidualZ / (1.5*ResidualAthrd+ra);
     209  ec = 1.44 * theZ * ResidualZ / (1.7*ResidualAthrd+ra);
     210  ecsq = ec * ec;
     211  p = p0 + p1/ec + p2/ecsq;
     212  landa = landa0*ResidualA + landa1;
     213  a = std::pow(ResidualA,mu1);
     214  mu = mu0 * a;
     215  nu = a* (nu0+nu1*ec+nu2*ecsq); 
     216  xnulam = nu / landa;
     217  if (xnulam > spill) xnulam=0.;
     218  if (xnulam >= flow) etest = 1.2 *std::sqrt(xnulam);
     219 
     220  a = -2.*p*ec + landa - nu/ecsq;
     221  b = p*ecsq + mu + 2.*nu/ec;
     222  ecut = 0.;
     223  cut = a*a - 4.*p*b;
     224  if (cut > 0.) ecut = std::sqrt(cut);
     225  ecut = (ecut-a) / (p+p);
     226  ecut2 = ecut;
     227  if (cut < 0.) ecut2 = ecut - 2.;
     228  elab = K * FragmentA / ResidualA;
     229  sig = 0.;
     230 
     231  if (elab <= ec) { //start for E<Ec
     232    if (elab > ecut2)  sig = (p*elab*elab+a*elab+b) * signor;
     233  }           //end for E<Ec
     234  else {           //start for E>Ec
     235    sig = (landa*elab+mu+nu/elab) * signor;
     236    geom = 0.;
     237    if (xnulam < flow || elab < etest) return sig;
     238    geom = std::sqrt(theA*K);
     239    geom = 1.23*ResidualAthrd + ra + 4.573/geom;
     240    geom = 31.416 * geom * geom;
     241    sig = std::max(geom,sig);
     242  }           //end for E>Ec
     243  return sig;
     244 
     245}
     246
     247//   ************************** end of cross sections *******************************

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

    r819 r962  
    2626//
    2727// $Id: G4DeuteronEvaporationChannel.cc,v 1.4 2006/06/29 20:10:21 gunter Exp $
    28 // GEANT4 tag $Name: $
     28// GEANT4 tag $Name: geant4-09-02-ref-02 $
    2929//
    3030// Hadronic Process: Nuclear De-excitations

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

    r819 r962  
    2424// ********************************************************************
    2525//
    26 //
    27 // $Id: G4DeuteronEvaporationProbability.cc,v 1.4 2006/06/29 20:10:23 gunter Exp $
    28 // GEANT4 tag $Name: geant4-09-01-patch-02 $
     26//J.M. Quesada (August2008). Based on:
    2927//
    3028// Hadronic Process: Nuclear De-excitations
    31 // by V. Lara (Nov 1999)
    32 //
     29// by V. Lara (Oct 1998)
     30//
     31// Modif (03 September 2008) by J. M. Quesada for external choice of inverse
     32// cross section option
    3333
    3434
    3535#include "G4DeuteronEvaporationProbability.hh"
    3636
     37
    3738G4DeuteronEvaporationProbability::G4DeuteronEvaporationProbability() :
    38     G4EvaporationProbability(2,1,6) // A,Z,Gamma
    39 {
    40     std::vector<G4double>::size_type NumExcitedStatesEnergy = 31+1;
    41     std::vector<G4int>::size_type NumExcitedStatesSpin = 31+1;
    42     ExcitEnergies.reserve(NumExcitedStatesEnergy);
    43     ExcitSpins.reserve(NumExcitedStatesSpin);
    44     ExcitEnergies.insert(ExcitEnergies.begin(),NumExcitedStatesEnergy,0.0);
    45     ExcitSpins.insert(ExcitSpins.begin(),NumExcitedStatesSpin,0);
    46 
    47 
    48     ExcitEnergies[15] = 6.18*MeV;
    49     ExcitEnergies[17] = 2.15*MeV;
    50     ExcitEnergies[18] = 5.02*MeV;
    51     ExcitEnergies[19] = 2.65*MeV;
    52     ExcitEnergies[20] = 4.80*MeV;
    53     ExcitEnergies[22] = 3.85*MeV;
    54     ExcitEnergies[23] = 6.96*MeV;
    55     ExcitEnergies[25] = 4.92*MeV;
    56     ExcitEnergies[26] = 7.22*MeV;
    57     ExcitEnergies[27] = 0.40*MeV;
    58     ExcitEnergies[28] = 6.83*MeV;
    59     ExcitEnergies[29] = 7.12*MeV;
    60     ExcitEnergies[30] = 3.84*MeV;
    61     ExcitEnergies[31] = 3.92*MeV;
    62 
    63     ExcitSpins[15] = 1;
    64     ExcitSpins[17] = 3;
    65     ExcitSpins[18] = 4;
    66     ExcitSpins[19] = 4;
    67     ExcitSpins[20] = 4;
    68     ExcitSpins[22] = 6;
    69     ExcitSpins[23] = 6;
    70     ExcitSpins[25] = 1;
    71     ExcitSpins[26] = 10;
    72     ExcitSpins[27] = 3;
    73     ExcitSpins[28] = 10;
    74     ExcitSpins[29] = 3;
    75     ExcitSpins[30] = 6;
    76     ExcitSpins[31] = 5;
    77        
    78     SetExcitationEnergiesPtr(&ExcitEnergies);
    79     SetExcitationSpinsPtr(&ExcitSpins);
    80        
    81 }
     39    G4EvaporationProbability(2,1,3,&theCoulombBarrier) // A,Z,Gamma (fixed JMQ)
     40{       }
    8241G4DeuteronEvaporationProbability::G4DeuteronEvaporationProbability(const G4DeuteronEvaporationProbability &) : G4EvaporationProbability()
    8342{
    8443    throw G4HadronicException(__FILE__, __LINE__, "G4DeuteronEvaporationProbability::copy_constructor meant to not be accessable");
    8544}
    86 
    87 
    8845
    8946
    … …  
    10158}
    10259
     60
    10361G4bool G4DeuteronEvaporationProbability::operator!=(const G4DeuteronEvaporationProbability &) const
    10462{
    … …  
    10664}
    10765
    108 G4double G4DeuteronEvaporationProbability::CCoeficient(const G4double aZ) const
     66
     67G4double G4DeuteronEvaporationProbability::CalcAlphaParam(const G4Fragment & fragment)
     68{
     69    return 1.0 + CCoeficient(static_cast<G4double>(fragment.GetZ()-GetZ()));
     70}
     71
     72
     73G4double G4DeuteronEvaporationProbability::CalcBetaParam(const G4Fragment & )
     74{
     75    return 0.0;
     76}
     77
     78
     79G4double G4DeuteronEvaporationProbability::CCoeficient(const G4double aZ)
    10980{
    11081    // Data comes from
    … …  
    12798       
    12899}
     100
     101
     102///////////////////////////////////////////////////////////////////////////////////
     103//J. M. Quesada (Dec 2007-June 2008): New inverse reaction cross sections
     104//OPT=0 Dostrovski's parameterization
     105//OPT=1,2 Chatterjee's paramaterization
     106//OPT=3,4 Kalbach's parameterization
     107//
     108G4double G4DeuteronEvaporationProbability::CrossSection(const  G4Fragment & fragment, const  G4double K)
     109{
     110       theA=GetA();
     111       theZ=GetZ();
     112       ResidualA=fragment.GetA()-theA;
     113       ResidualZ=fragment.GetZ()-theZ;
     114 
     115       ResidualAthrd=std::pow(ResidualA,0.33333);
     116       FragmentA=fragment.GetA();
     117       FragmentAthrd=std::pow(FragmentA,0.33333);
     118
     119       if (OPTxs==0) {std::ostringstream errOs;
     120         errOs << "We should'n be here (OPT =0) at evaporation cross section calculation (deuterons)!!"  <<G4endl;
     121         throw G4HadronicException(__FILE__, __LINE__, errOs.str());
     122         return 0.;}
     123       if( OPTxs==1 || OPTxs==2) return G4DeuteronEvaporationProbability::GetOpt12( K);
     124       else if (OPTxs==3 || OPTxs==4)  return G4DeuteronEvaporationProbability::GetOpt34( K);
     125       else{
     126         std::ostringstream errOs;
     127         errOs << "BAD Deuteron CROSS SECTION OPTION AT EVAPORATION!!"  <<G4endl;
     128         throw G4HadronicException(__FILE__, __LINE__, errOs.str());
     129         return 0.;
     130       }
     131}
     132
     133
     134//********************* OPT=1,2 : Chatterjee's cross section ************************
     135//(fitting to cross section from Bechetti & Greenles OM potential)
     136
     137G4double G4DeuteronEvaporationProbability::GetOpt12(const  G4double K)
     138{
     139
     140  G4double Kc=K;
     141
     142// JMQ xsec is set constat above limit of validity
     143 if (K>50) Kc=50;
     144
     145 G4double landa ,mu ,nu ,p , Ec,q,r,ji,xs;
     146
     147 
     148 G4double    p0 = -38.21;
     149 G4double    p1 = 922.6;
     150 G4double    p2 = -2804.;
     151 G4double    landa0 = -0.0323;
     152 G4double    landa1 = -5.48;
     153 G4double    mu0 = 336.1;
     154 G4double    mu1 = 0.48;
     155 G4double    nu0 = 524.3;
     156 G4double    nu1 = -371.8;
     157 G4double    nu2 = -5.924; 
     158 G4double    delta=1.2;           
     159 
     160
     161 Ec = 1.44*theZ*ResidualZ/(1.5*ResidualAthrd+delta);
     162 p = p0 + p1/Ec + p2/(Ec*Ec);
     163 landa = landa0*ResidualA + landa1;
     164 mu = mu0*std::pow(ResidualA,mu1);
     165 nu = std::pow(ResidualA,mu1)*(nu0 + nu1*Ec + nu2*(Ec*Ec));
     166 q = landa - nu/(Ec*Ec) - 2*p*Ec;
     167 r = mu + 2*nu/Ec + p*(Ec*Ec);
     168 
     169 ji=std::max(Kc,Ec);
     170 if(Kc < Ec) { xs = p*Kc*Kc + q*Kc + r;}
     171 else {xs = p*(Kc - ji)*(Kc - ji) + landa*Kc + mu + nu*(2 - Kc/ji)/ji ;}
     172                 
     173 if (xs <0.0) {xs=0.0;}
     174 
     175 return xs;
     176
     177}
     178
     179
     180// *********** OPT=3,4 : Kalbach's cross sections (from PRECO code)*************
     181G4double G4DeuteronEvaporationProbability::GetOpt34(const  G4double K)
     182//     ** d from o.m. of perey and perey
     183{
     184
     185  G4double landa, mu, nu, p , signor(1.),sig;
     186  G4double ec,ecsq,xnulam,etest(0.),a;
     187  G4double b,ecut,cut,ecut2,geom,elab;
     188
     189
     190  G4double     flow = 1.e-18;
     191  G4double     spill= 1.e+18;
     192
     193
     194  G4double     p0 = 0.798;
     195  G4double     p1 = 420.3;
     196  G4double     p2 = -1651.;
     197  G4double     landa0 = 0.00619;
     198  G4double     landa1 = -7.54;
     199  G4double     mu0 = 583.5;
     200  G4double     mu1 = 0.337;
     201  G4double     nu0 = 421.8;
     202  G4double     nu1 = -474.5;
     203  G4double     nu2 = -3.592;     
     204 
     205  G4double      ra=0.80;
     206       
     207  //JMQ 13/02/09 increase of reduced radius to lower the barrier
     208  // ec = 1.44 * theZ * ResidualZ / (1.5*ResidualAthrd+ra);
     209  ec = 1.44 * theZ * ResidualZ / (1.7*ResidualAthrd+ra);
     210  ecsq = ec * ec;
     211  p = p0 + p1/ec + p2/ecsq;
     212  landa = landa0*ResidualA + landa1;
     213  a = std::pow(ResidualA,mu1);
     214  mu = mu0 * a;
     215  nu = a* (nu0+nu1*ec+nu2*ecsq); 
     216  xnulam = nu / landa;
     217  if (xnulam > spill) xnulam=0.;
     218  if (xnulam >= flow) etest = 1.2 *std::sqrt(xnulam);
     219
     220  a = -2.*p*ec + landa - nu/ecsq;
     221  b = p*ecsq + mu + 2.*nu/ec;
     222  ecut = 0.;
     223  cut = a*a - 4.*p*b;
     224  if (cut > 0.) ecut = std::sqrt(cut);
     225  ecut = (ecut-a) / (p+p);
     226  ecut2 = ecut;
     227  if (cut < 0.) ecut2 = ecut - 2.;
     228  elab = K * FragmentA / ResidualA;
     229  sig = 0.;
     230
     231  if (elab <= ec) { //start for E<Ec
     232    if (elab > ecut2)  sig = (p*elab*elab+a*elab+b) * signor;
     233  }           //end for E<Ec
     234  else {           //start for E>Ec
     235    sig = (landa*elab+mu+nu/elab) * signor;
     236    geom = 0.;
     237    if (xnulam < flow || elab < etest) return sig;
     238    geom = std::sqrt(theA*K);
     239    geom = 1.23*ResidualAthrd + ra + 4.573/geom;
     240    geom = 31.416 * geom * geom;
     241    sig = std::max(geom,sig);
     242  }           //end for E>Ec
     243  return sig;
     244 
     245}
     246
     247//   ************************** end of cross sections *******************************

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

    r819 r962  
    2525//
    2626//
    27 // $Id: G4Evaporation.cc,v 1.7 2007/02/14 13:37:49 ahoward Exp $
    28 // GEANT4 tag $Name: $
     27// $Id: G4Evaporation.cc,v 1.12 2008/12/09 17:57:36 ahoward Exp $
     28// GEANT4 tag $Name: geant4-09-02-ref-02 $
    2929//
    3030// Hadronic Process: Nuclear De-excitations
    … …  
    3333// Alex Howard - added protection for negative probabilities in the sum, 14/2/07
    3434//
     35// Modif (03 September 2008) by J. M. Quesada for external choice of inverse
     36// cross section option
     37// JMQ (06 September 2008) Also external choices have been added for
     38// superimposed Coulomb barrier (if useSICBis set true, by default is false)
     39
    3540#include "G4Evaporation.hh"
    3641#include "G4EvaporationFactory.hh"
    … …  
    8994}
    9095
    91 
    9296G4FragmentVector * G4Evaporation::BreakItUp(const G4Fragment &theNucleus)
    9397{
    … …  
    105109    // Number of channels
    106110    G4int TotNumberOfChannels = theChannels->size(); 
    107        
    108111
    109112    // Starts loop over evaporated particles
    110113    for (;;)
    111       {
     114 
     115   {   
    112116        // loop over evaporation channels
    113117        std::vector<G4VEvaporationChannel*>::iterator i;
    114118        for (i=theChannels->begin(); i != theChannels->end(); i++)
    115119          {
     120  // for inverse cross section choice
     121            (*i)->SetOPTxs(OPTxs);
     122  // for superimposed Coulomb Barrier for inverse cross sections
     123            (*i)->UseSICB(useSICB);
     124
    116125            (*i)->Initialize(theResidualNucleus);
    117126          }
    … …  
    172181            if( j >= TotNumberOfChannels )
    173182              {
     183                G4cerr << " Residual A: " << theResidualNucleus.GetA() << " Residual Z: " << theResidualNucleus.GetZ() << " Excitation Energy: " << theResidualNucleus.GetExcitationEnergy() << G4endl;
     184                G4cerr << " j has not chosen a channel, j = " << j << " TotNumberOfChannels " << TotNumberOfChannels << " Total Probability: " << TotalProbability << G4endl;
     185                for (j=0; j < TotNumberOfChannels; j++)
     186                  {
     187                    G4cerr << " j: " << j << " EmissionProbChannel: " << EmissionProbChannel[j] << " and shoot: " << shoot << " (<ProbChannel?) " << G4endl;
     188                  }             
    174189                throw G4HadronicException(__FILE__, __LINE__,  "G4Evaporation::BreakItUp: Can't define emission probability of the channels!" );
    175190              }
    … …  
    216231                  theResidualNucleus.SetCreatorModel(G4String("ResidualNucleus"));
    217232#endif
     233
    218234                }
    219235              }

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

    r819 r962  
    2525//
    2626//
    27 // $Id: G4EvaporationChannel.cc,v 1.5 2006/06/29 20:10:27 gunter Exp $
    28 // GEANT4 tag $Name:  $
     27//J.M. Quesada (August2008). Based on:
    2928//
    3029// Hadronic Process: Nuclear De-excitations
    3130// by V. Lara (Oct 1998)
    3231//
     32// Modif (03 September 2008) by J. M. Quesada for external choice of inverse
     33// cross section option
     34// JMQ (06 September 2008) Also external choices have been added for
     35// superimposed Coulomb barrier (if useSICB is set true, by default is false)
     36
    3337
    3438#include "G4EvaporationChannel.hh"
    … …  
    3640
    3741
    38 G4EvaporationChannel::G4EvaporationChannel(const G4int theA, const G4int theZ,
     42
     43G4EvaporationChannel::G4EvaporationChannel(const G4int anA, const G4int aZ, const G4String & aName,
    3944                                           G4VEmissionProbability * aEmissionStrategy,
    40                                            G4VCoulombBarrier * aCoulombBarrier):
    41     A(theA),
    42     Z(theZ),
     45                                           G4VCoulombBarrier * aCoulombBarrier):
     46    G4VEvaporationChannel(aName),
     47    theA(anA),
     48    theZ(aZ),
    4349    theEvaporationProbabilityPtr(aEmissionStrategy),
    4450    theCoulombBarrierPtr(aCoulombBarrier),
    … …  
    4753{
    4854    theLevelDensityPtr = new G4EvaporationLevelDensityParameter;
    49     MyOwnLevelDensity = true;
    50 }
    51 
    52 G4EvaporationChannel::G4EvaporationChannel(const G4int theA, const G4int theZ, const G4String & aName,
    53                                            G4VEmissionProbability * aEmissionStrategy,
    54                                            G4VCoulombBarrier * aCoulombBarrier):
    55     G4VEvaporationChannel(aName),
    56     A(theA),
    57     Z(theZ),
    58     theEvaporationProbabilityPtr(aEmissionStrategy),
    59     theCoulombBarrierPtr(aCoulombBarrier),
    60     EmissionProbability(0.0),
    61     MaximalKineticEnergy(-1000.0)
    62 {
    63     theLevelDensityPtr = new G4EvaporationLevelDensityParameter;
    64     MyOwnLevelDensity = true;
    65 }
    66 
    67 G4EvaporationChannel::G4EvaporationChannel(const G4int theA, const G4int theZ, const G4String * aName,
    68                                            G4VEmissionProbability * aEmissionStrategy,
    69                                            G4VCoulombBarrier * aCoulombBarrier):
    70     G4VEvaporationChannel(aName),
    71     A(theA),
    72     Z(theZ),
    73     theEvaporationProbabilityPtr(aEmissionStrategy),
    74     theCoulombBarrierPtr(aCoulombBarrier),
    75     EmissionProbability(0.0),
    76     MaximalKineticEnergy(-1000.0)
    77 {
    78     theLevelDensityPtr = new G4EvaporationLevelDensityParameter;
    79     MyOwnLevelDensity = true;
    80 }
    81 
     55//    MyOwnLevelDensity = true; 
     56}
    8257
    8358G4EvaporationChannel::~G4EvaporationChannel()
    8459{
    85     if (MyOwnLevelDensity) delete theLevelDensityPtr;
    86 }
    87 
    88 
    89 
     60  delete theLevelDensityPtr;
     61}
    9062
    9163G4EvaporationChannel::G4EvaporationChannel(const G4EvaporationChannel & ) : G4VEvaporationChannel()
    … …  
    11284}
    11385
    114 
     86//void G4EvaporationChannel::SetLevelDensityParameter(G4VLevelDensityParameter * aLevelDensity)
     87//  {
     88//    if (MyOwnLevelDensity) delete theLevelDensityPtr;
     89//    theLevelDensityPtr = aLevelDensity;
     90//    MyOwnLevelDensity = false;
     91//  }
    11592
    11693void G4EvaporationChannel::Initialize(const G4Fragment & fragment)
    11794{
    118 
    119     G4int anA = static_cast<G4int>(fragment.GetA());
    120     G4int aZ = static_cast<G4int>(fragment.GetZ());
    121     AResidual = anA - A;
    122     ZResidual = aZ - Z;
    123 
    124     // Effective excitation energy
    125     G4double ExEnergy = fragment.GetExcitationEnergy() -
    126       G4PairingCorrection::GetInstance()->GetPairingCorrection(anA,aZ);
    127 
    128     // We only take into account channels which are physically allowed
    129     if (AResidual <= 0 || ZResidual <= 0 || AResidual < ZResidual ||
    130         (AResidual == ZResidual && AResidual > 1) || ExEnergy <= 0.0) {
    131         //              LevelDensityParameter = 0.0;
    132         CoulombBarrier = 0.0;
    133         //              BindingEnergy = 0.0;
    134         MaximalKineticEnergy = -1000.0*MeV;
    135         EmissionProbability = 0.0;
    136     } else {
    137         // // Get Level Density
    138         // LevelDensityParameter = theLevelDensityPtr->LevelDensityParameter(anA,aZ,ExEnergy);
    139 
    140         // Coulomb Barrier calculation
    141         CoulombBarrier = theCoulombBarrierPtr->GetCoulombBarrier(AResidual,ZResidual,ExEnergy);
    142        
    143         // // Binding Enegy (for separate fragment from nucleus)
    144         // BindingEnergy = CalcBindingEnergy(anA,aZ);
    145 
    146         // Maximal Kinetic Energy
    147         MaximalKineticEnergy = CalcMaximalKineticEnergy(G4ParticleTable::GetParticleTable()->
    148                                                         GetIonTable()->GetNucleusMass(aZ,anA)+ExEnergy);
    149                
    150         // Emission probability
    151         if (MaximalKineticEnergy <= 0.0) EmissionProbability = 0.0;
    152         else {
    153             // Total emission probability for this channel
    154             EmissionProbability = theEvaporationProbabilityPtr->EmissionProbability(fragment,MaximalKineticEnergy);
    155 
    156         }
     95  //for inverse cross section choice
     96  theEvaporationProbabilityPtr->SetOPTxs(OPTxs);
     97  // for superimposed Coulomb Barrier for inverse cross sections
     98  theEvaporationProbabilityPtr->UseSICB(useSICB);
     99 
     100 
     101  G4int FragmentA = static_cast<G4int>(fragment.GetA());
     102  G4int FragmentZ = static_cast<G4int>(fragment.GetZ());
     103  ResidualA = FragmentA - theA;
     104  ResidualZ = FragmentZ - theZ;
     105 
     106  //Effective excitation energy
     107  G4double ExEnergy = fragment.GetExcitationEnergy() -
     108    G4PairingCorrection::GetInstance()->GetPairingCorrection(FragmentA,FragmentZ);
     109 
     110 
     111  // Only channels which are physically allowed are taken into account
     112  if (ResidualA <= 0 || ResidualZ <= 0 || ResidualA < ResidualZ ||
     113      (ResidualA == ResidualZ && ResidualA > 1) || ExEnergy <= 0.0) {
     114    CoulombBarrier=0.0;
     115    MaximalKineticEnergy = -1000.0*MeV;
     116    EmissionProbability = 0.0;
     117  } else {
     118    CoulombBarrier = theCoulombBarrierPtr->GetCoulombBarrier(ResidualA,ResidualZ,ExEnergy);
     119    // Maximal Kinetic Energy
     120    MaximalKineticEnergy = CalcMaximalKineticEnergy
     121      (G4ParticleTable::GetParticleTable()->
     122       GetIonTable()->GetNucleusMass(FragmentZ,FragmentA)+ExEnergy);
     123   
     124    // Emission probability
     125    // Protection for the case Tmax<V. If not set in this way we could end up in an
     126    // infinite loop in  the method GetKineticEnergy if OPTxs!=0 && useSICB=true.
     127    // Of course for OPTxs=0 we have the Coulomb barrier
     128   
     129    G4double limit;
     130    if (OPTxs==0 || (OPTxs!=0 && useSICB))
     131      limit= CoulombBarrier;
     132    else limit=0.;
     133 
     134    // The threshold for charged particle emission must be  set to 0 if Coulomb
     135    //cutoff  is included in the cross sections
     136    //if (MaximalKineticEnergy <= 0.0) EmissionProbability = 0.0; 
     137    if (MaximalKineticEnergy <= limit) EmissionProbability = 0.0;
     138    else {
     139      // Total emission probability for this channel
     140      EmissionProbability = theEvaporationProbabilityPtr->
     141        EmissionProbability(fragment, MaximalKineticEnergy);
    157142    }
    158        
    159     return;
    160 }
    161 
     143  }
     144 
     145  return;
     146}
    162147
    163148G4FragmentVector * G4EvaporationChannel::BreakUp(const G4Fragment & theNucleus)
    164149{
    165 
    166     G4double EvaporatedKineticEnergy = CalcKineticEnergy();
    167     G4double EvaporatedMass = G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(Z,A);
    168     G4double EvaporatedEnergy = EvaporatedKineticEnergy + EvaporatedMass;
    169 
    170     G4ThreeVector momentum(IsotropicVector(std::sqrt(EvaporatedKineticEnergy*
    171                                                 (EvaporatedKineticEnergy+2.0*EvaporatedMass))));
    172  
    173     momentum.rotateUz(theNucleus.GetMomentum().vect().unit());
    174 
    175     G4LorentzVector EvaporatedMomentum(momentum,EvaporatedEnergy);
    176     EvaporatedMomentum.boost(theNucleus.GetMomentum().boostVector());
    177 
    178     G4Fragment * EvaporatedFragment = new G4Fragment(A,Z,EvaporatedMomentum);
     150  G4double EvaporatedKineticEnergy=GetKineticEnergy(theNucleus);
     151 
     152  G4double EvaporatedMass = G4ParticleTable::GetParticleTable()->GetIonTable()->
     153    GetIonMass(theZ,theA);
     154  G4double EvaporatedEnergy = EvaporatedKineticEnergy + EvaporatedMass;
     155 
     156  G4ThreeVector momentum(IsotropicVector
     157                         (std::sqrt(EvaporatedKineticEnergy*
     158                                    (EvaporatedKineticEnergy+2.0*EvaporatedMass))));
     159 
     160  momentum.rotateUz(theNucleus.GetMomentum().vect().unit());
     161 
     162  G4LorentzVector EvaporatedMomentum(momentum,EvaporatedEnergy);
     163  EvaporatedMomentum.boost(theNucleus.GetMomentum().boostVector());
     164 
     165  G4Fragment * EvaporatedFragment = new G4Fragment(theA,theZ,EvaporatedMomentum);
    179166#ifdef PRECOMPOUND_TEST
    180     EvaporatedFragment->SetCreatorModel(G4String("G4Evaporation"));
     167  EvaporatedFragment->SetCreatorModel(G4String("G4Evaporation"));
    181168#endif
    182     // ** And now the residual nucleus **
    183     G4double theExEnergy = theNucleus.GetExcitationEnergy();
    184     G4double theMass = G4ParticleTable::GetParticleTable()->GetIonTable()->
    185         GetNucleusMass(static_cast<G4int>(theNucleus.GetZ()),
    186                        static_cast<G4int>(theNucleus.GetA()));
    187     G4double ResidualEnergy = theMass + (theExEnergy - EvaporatedKineticEnergy) - EvaporatedMass;
    188        
    189     G4LorentzVector ResidualMomentum(-momentum,ResidualEnergy);
    190     ResidualMomentum.boost(theNucleus.GetMomentum().boostVector());
    191        
    192     G4Fragment * ResidualFragment = new G4Fragment( AResidual, ZResidual, ResidualMomentum );
     169  // ** And now the residual nucleus **
     170  G4double theExEnergy = theNucleus.GetExcitationEnergy();
     171  G4double theMass = G4ParticleTable::GetParticleTable()->GetIonTable()->
     172    GetNucleusMass(static_cast<G4int>(theNucleus.GetZ()),
     173                   static_cast<G4int>(theNucleus.GetA()));
     174  G4double ResidualEnergy = theMass +
     175    (theExEnergy - EvaporatedKineticEnergy) - EvaporatedMass;
     176 
     177  G4LorentzVector ResidualMomentum(-momentum,ResidualEnergy);
     178  ResidualMomentum.boost(theNucleus.GetMomentum().boostVector());
     179 
     180  G4Fragment * ResidualFragment = new G4Fragment(ResidualA, ResidualZ, ResidualMomentum );
     181 
    193182#ifdef PRECOMPOUND_TEST
    194     ResidualFragment->SetCreatorModel(G4String("ResidualNucleus"));
     183  ResidualFragment->SetCreatorModel(G4String("ResidualNucleus"));
    195184#endif
    196     G4FragmentVector * theResult = new G4FragmentVector;
    197 
     185  G4FragmentVector * theResult = new G4FragmentVector;
     186 
    198187#ifdef debug
    199     G4double Efinal = ResidualMomentum.e() + EvaporatedMomentum.e();
    200     G4ThreeVector Pfinal = ResidualMomentum.vect() + EvaporatedMomentum.vect();
    201     if (std::abs(Efinal-theNucleus.GetMomentum().e()) > 1.0*keV) {
    202         G4cout << "@@@@@@@@@@@@@@@@@@@@@ G4Evaporation Chanel: ENERGY @@@@@@@@@@@@@@@@" << G4endl;
    203         G4cout << "Initial : " << theNucleus.GetMomentum().e()/MeV << " MeV    Final :"
    204                <<Efinal/MeV << " MeV    Delta: " <<  (Efinal-theNucleus.GetMomentum().e())/MeV
    205                << " MeV" << G4endl;
    206     }
    207     if (std::abs(Pfinal.x()-theNucleus.GetMomentum().x()) > 1.0*keV ||
    208         std::abs(Pfinal.y()-theNucleus.GetMomentum().y()) > 1.0*keV ||
    209         std::abs(Pfinal.z()-theNucleus.GetMomentum().z()) > 1.0*keV ) {
    210         G4cout << "@@@@@@@@@@@@@@@@@@@@@ G4Evaporation Chanel: MOMENTUM @@@@@@@@@@@@@@@@" << G4endl;
    211         G4cout << "Initial : " << theNucleus.GetMomentum().vect() << " MeV    Final :"
    212                <<Pfinal/MeV << " MeV    Delta: " <<  Pfinal-theNucleus.GetMomentum().vect()
    213                << " MeV" << G4endl;   
    214 
    215     }
     188  G4double Efinal = ResidualMomentum.e() + EvaporatedMomentum.e();
     189  G4ThreeVector Pfinal = ResidualMomentum.vect() + EvaporatedMomentum.vect();
     190  if (std::abs(Efinal-theNucleus.GetMomentum().e()) > 1.0*keV) {
     191    G4cout << "@@@@@@@@@@@@@@@@@@@@@ G4Evaporation Chanel: ENERGY @@@@@@@@@@@@@@@@" << G4endl;
     192    G4cout << "Initial : " << theNucleus.GetMomentum().e()/MeV << " MeV    Final :"
     193           <<Efinal/MeV << " MeV    Delta: " <<  (Efinal-theNucleus.GetMomentum().e())/MeV
     194           << " MeV" << G4endl;
     195  }
     196  if (std::abs(Pfinal.x()-theNucleus.GetMomentum().x()) > 1.0*keV ||
     197      std::abs(Pfinal.y()-theNucleus.GetMomentum().y()) > 1.0*keV ||
     198      std::abs(Pfinal.z()-theNucleus.GetMomentum().z()) > 1.0*keV ) {
     199    G4cout << "@@@@@@@@@@@@@@@@@@@@@ G4Evaporation Chanel: MOMENTUM @@@@@@@@@@@@@@@@" << G4endl;
     200    G4cout << "Initial : " << theNucleus.GetMomentum().vect() << " MeV    Final :"
     201           <<Pfinal/MeV << " MeV    Delta: " <<  Pfinal-theNucleus.GetMomentum().vect()
     202           << " MeV" << G4endl;   
     203   
     204  }
    216205#endif
    217     theResult->push_back(EvaporatedFragment);
    218     theResult->push_back(ResidualFragment);
    219     return theResult;
     206  theResult->push_back(EvaporatedFragment);
     207  theResult->push_back(ResidualFragment);
     208  return theResult;
    220209}
    221210
    222 
    223 
    224 // G4double G4EvaporationChannel::CalcBindingEnergy(const G4int anA, const G4int aZ)
    225 // // Calculate Binding Energy for separate fragment from nucleus
    226 // {
    227 //      // Mass Excess for residual nucleus
    228 //      G4double ResNucMassExcess = G4NucleiProperties::GetNuclearMass(AResidual,ZResidual);
    229 //      // Mass Excess for fragment
    230 //      G4double FragmentMassExcess = G4NucleiProperties::GetNuclearMass(A,Z);
    231 //      // Mass Excess for Compound Nucleus
    232 //      G4double NucleusMassExcess = G4NucleiProperties::GetNuclearMass(anA,aZ);
    233 //
    234 //      return ResNucMassExcess + FragmentMassExcess - NucleusMassExcess;
    235 // }
    236 
    237 
     211/////////////////////////////////////////
     212// Calculates the maximal kinetic energy that can be carried by fragment.
    238213G4double G4EvaporationChannel::CalcMaximalKineticEnergy(const G4double NucleusTotalE)
    239     // Calculate maximal kinetic energy that can be carried by fragment.
    240 {
    241     G4double ResidualMass = G4ParticleTable::GetParticleTable()->GetIonTable()->GetNucleusMass( ZResidual, AResidual );
    242     G4double EvaporatedMass = G4ParticleTable::GetParticleTable()->GetIonTable()->GetNucleusMass( Z, A );
    243        
    244     G4double T = (NucleusTotalE*NucleusTotalE + EvaporatedMass*EvaporatedMass - ResidualMass*ResidualMass)/
    245         (2.0*NucleusTotalE) -
    246         EvaporatedMass - CoulombBarrier;
    247        
    248     return T;
    249 }
    250 
    251 
    252 
    253 
    254 G4double G4EvaporationChannel::CalcKineticEnergy(void)
    255     // Samples fragment kinetic energy.
     214{
     215  G4double ResidualMass = G4ParticleTable::GetParticleTable()->
     216    GetIonTable()->GetNucleusMass( ResidualZ, ResidualA );
     217  G4double EvaporatedMass = G4ParticleTable::GetParticleTable()->
     218    GetIonTable()->GetNucleusMass( theZ, theA );
     219
     220  // This is the "true" assimptotic kinetic energy (from energy conservation)   
     221  G4double Tmax = (NucleusTotalE*NucleusTotalE + EvaporatedMass*EvaporatedMass -               
     222                   ResidualMass*ResidualMass)/(2.0*NucleusTotalE) - EvaporatedMass;
     223 
     224  //JMQ (13-09-08) bug fixed: in the original version the Tmax is calculated
     225  //at the Coulomb barrier
     226  //IMPORTANT: meaning of Tmax differs in OPTxs=0 and OPTxs!=0
     227  //When OPTxs!=0 Tmax is the TRUE (assimptotic) maximal kinetic energy
     228 
     229  if(OPTxs==0)
     230    Tmax=Tmax- CoulombBarrier;
     231 
     232  return Tmax;
     233}
     234
     235///////////////////////////////////////////
     236//JMQ: New method for MC sampling of kinetic energy. Substitutes old CalcKineticEnergy
     237G4double G4EvaporationChannel::GetKineticEnergy(const G4Fragment & aFragment)
     238{
     239 
     240  if (OPTxs==0) {
    256241    // It uses Dostrovsky's approximation for the inverse reaction cross
    257     // in the probability for fragment emisson
    258 {
    259     if (MaximalKineticEnergy < 0.0)
    260         throw G4HadronicException(__FILE__, __LINE__, "G4EvaporationChannel::CalcKineticEnergy: maximal kinetic energy is less than 0");
    261 
    262     G4double Rb = 4.0*theLevelDensityPtr->LevelDensityParameter(AResidual+A,ZResidual+Z,MaximalKineticEnergy)*
    263         MaximalKineticEnergy;
     242    // in the probability for fragment emission
     243    // MaximalKineticEnergy energy in the original version (V.Lara) was calculated at
     244    //the Coulomb barrier.
     245   
     246   
     247    if (MaximalKineticEnergy < 0.0)   
     248      throw G4HadronicException(__FILE__, __LINE__,
     249                                "G4EvaporationChannel::CalcKineticEnergy: maximal kinetic at the Coulomb barrier is less than 0");
     250   
     251    G4double Rb = 4.0*theLevelDensityPtr->
     252      LevelDensityParameter(ResidualA+theA,ResidualZ+theZ,MaximalKineticEnergy)*
     253      MaximalKineticEnergy;
    264254    G4double RbSqrt = std::sqrt(Rb);
    265255    G4double PEX1 = 0.0;
    … …  
    268258    G4double FRk = 0.0;
    269259    do {
    270         G4double RandNumber = G4UniformRand();
    271         Rk = 1.0 + (1./RbSqrt)*std::log(RandNumber + (1.0-RandNumber)*PEX1);
    272         G4double Q1 = 1.0;
    273         G4double Q2 = 1.0;
    274         if (Z == 0) { // for emitted neutron
    275             G4double Beta = (2.12/std::pow(G4double(AResidual),2./3.) - 0.05)*MeV/
    276                 (0.76 + 2.2/std::pow(G4double(AResidual),1./3.));
    277             Q1 = 1.0 + Beta/(MaximalKineticEnergy);
    278             Q2 = Q1*std::sqrt(Q1);
    279         }
    280    
    281         FRk = (3.0*std::sqrt(3.0)/2.0)/Q2 * Rk * (Q1 - Rk*Rk);
    282    
     260      G4double RandNumber = G4UniformRand();
     261      Rk = 1.0 + (1./RbSqrt)*std::log(RandNumber + (1.0-RandNumber)*PEX1);
     262      G4double Q1 = 1.0;
     263      G4double Q2 = 1.0;
     264      if (theZ == 0) { // for emitted neutron
     265        G4double Beta = (2.12/std::pow(G4double(ResidualA),2./3.) - 0.05)*MeV/
     266          (0.76 + 2.2/std::pow(G4double(ResidualA),1./3.));
     267        Q1 = 1.0 + Beta/(MaximalKineticEnergy);
     268        Q2 = Q1*std::sqrt(Q1);
     269      }
     270     
     271      FRk = (3.0*std::sqrt(3.0)/2.0)/Q2 * Rk * (Q1 - Rk*Rk);
     272     
    283273    } while (FRk < G4UniformRand());
    284 
     274   
    285275    G4double result =  MaximalKineticEnergy * (1.0-Rk*Rk) + CoulombBarrier;
    286 
    287276    return result;
    288 }
    289  
     277  } else if (OPTxs==1 || OPTxs==2 || OPTxs==3 || OPTxs==4) {
     278   
     279    G4double V;
     280    if(useSICB) V= CoulombBarrier;
     281    else V=0.;
     282    //If Coulomb barrier is just included  in the cross sections
     283    //  G4double V=0.;
     284
     285    G4double Tmax=MaximalKineticEnergy;
     286    G4double T(0.0);
     287    G4double NormalizedProbability(1.0);
     288   
     289    // A pointer is created in order to access the distribution function.
     290    G4EvaporationProbability * G4EPtemp;
     291   
     292    if (theA==1 && theZ==0) G4EPtemp=new G4NeutronEvaporationProbability();
     293    else if (theA==1 && theZ==1) G4EPtemp=new G4ProtonEvaporationProbability();
     294    else if (theA==2 && theZ==1 ) G4EPtemp=new G4DeuteronEvaporationProbability();
     295    else if (theA==3 && theZ==1 ) G4EPtemp=new G4TritonEvaporationProbability();
     296    else if (theA==3 && theZ==2 ) G4EPtemp=new G4He3EvaporationProbability();
     297    else if (theA==4 && theZ==2) G4EPtemp=new G4AlphaEvaporationProbability();
     298    else {
     299      std::ostringstream errOs;
     300      errOs << "ejected particle out of range in G4EvaporationChannel"  << G4endl;
     301      throw G4HadronicException(__FILE__, __LINE__, errOs.str());
     302    }
     303
     304      //for cross section selection and superimposed Coulom Barrier for xs
     305      G4EPtemp->SetOPTxs(OPTxs);
     306      G4EPtemp->UseSICB(useSICB);
     307
     308    do
     309      { 
     310        T=V+G4UniformRand()*(Tmax-V);
     311        NormalizedProbability=G4EPtemp->ProbabilityDistributionFunction(aFragment,T)/
     312          (this->GetEmissionProbability());
     313       
     314      }
     315    while (G4UniformRand() > NormalizedProbability);
     316   delete G4EPtemp;
     317   return T;
     318  } else{
     319    std::ostringstream errOs;
     320    errOs << "Bad option for energy sampling in evaporation"  <<G4endl;
     321    throw G4HadronicException(__FILE__, __LINE__, errOs.str());
     322  }
     323}
    290324
    291325G4ThreeVector G4EvaporationChannel::IsotropicVector(const G4double Magnitude)
    … …  
    300334                         Magnitude*CosTheta);
    301335    return Vector;
    302 }
    303 
    304 
    305 
     336            }
     337
     338
     339
     340
     341
     342   
     343
     344
     345 
     346

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

    r819 r962  
    2626//
    2727// $Id: G4EvaporationFactory.cc,v 1.4 2006/06/29 20:10:29 gunter Exp $
    28 // GEANT4 tag $Name: $
     28// GEANT4 tag $Name: geant4-09-02-ref-02 $
    2929//
    3030// Hadronic Process: Nuclear De-excitations

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

    r819 r962  
    2424// ********************************************************************
    2525//
    26 //
    27 // $Id: G4EvaporationProbability.cc,v 1.5 2006/06/29 20:10:31 gunter Exp $
    28 // GEANT4 tag $Name: geant4-09-01-patch-02 $
     26//J.M. Quesada (August2008). Based on:
    2927//
    3028// Hadronic Process: Nuclear De-excitations
    3129// by V. Lara (Oct 1998)
    3230//
    33 
     31// Modif (03 September 2008) by J. M. Quesada for external choice of inverse
     32// cross section option
     33// JMQ (06 September 2008) Also external choices have been added for
     34// superimposed Coulomb barrier (if useSICB is set true, by default is false)
     35//
     36// JMQ (14 february 2009) bug fixed in emission width: hbarc instead of hbar_Planck in the denominator
     37//
     38#include <iostream>
     39using namespace std;
    3440
    3541#include "G4EvaporationProbability.hh"
    … …  
    6369    return true;
    6470}
    65 
     71 
    6672G4double G4EvaporationProbability::EmissionProbability(const G4Fragment & fragment, const G4double anEnergy)
    6773{
    6874    G4double EmissionProbability = 0.0;
    6975    G4double MaximalKineticEnergy = anEnergy;
     76
    7077    if (MaximalKineticEnergy > 0.0 && fragment.GetExcitationEnergy() > 0.0) {
    71         EmissionProbability = CalcProbability(fragment,MaximalKineticEnergy);
    72         //              // Next there is a loop over excited states for this channel summing probabilities
    73         //              G4double SavedGamma = Gamma;
    74         //              for (G4int i = 0; i < ExcitationEnergies->length(); i++) {
    75         //                      if (ExcitationSpins->operator()(i) < 0.1) continue;
    76         //                      Gamma = ExcitationSpins->operator()(i)*A; // A is the channel mass number
    77         //                      // substract excitation energies
    78         //                      MaximalKineticEnergy -= ExcitationEnergies->operator()(i);
    79         //                      // update probability
    80         //                      EmissionProbability += CalcProbability(fragment,MaximalKineticEnergy);
    81         //                              EmissionProbability += tmp;
    82         //                      }
    83         //              // restore Gamma and MaximalKineticEnergy
    84         //              MaximalKineticEnergy = SavedMaximalKineticEnergy;
    85         //              Gamma = SavedGamma;
    86         //              }
     78        EmissionProbability = CalculateProbability(fragment, MaximalKineticEnergy);
     79
    8780    }
    8881    return EmissionProbability;
    8982}
    9083
    91 G4double G4EvaporationProbability::CalcProbability(const G4Fragment & fragment,
    92                                                    const G4double MaximalKineticEnergy)
    93     // Calculate integrated probability (width) for rvaporation channel
    94 {       
    95     G4double ResidualA = static_cast<G4double>(fragment.GetA() - theA);
    96     G4double ResidualZ = static_cast<G4double>(fragment.GetZ() - theZ);
     84////////////////////////////////////
     85
     86// Computes the integrated probability of evaporation channel
     87G4double G4EvaporationProbability::CalculateProbability(const G4Fragment & fragment, const G4double MaximalKineticEnergy)
     88{
     89    G4double ResidualA = fragment.GetA() - theA;
     90    G4double ResidualZ = fragment.GetZ() - theZ;
    9791    G4double U = fragment.GetExcitationEnergy();
     92   
     93 if (OPTxs==0) {
     94
    9895       
    9996    G4double NuclearMass = G4ParticleTable::GetParticleTable()->GetIonTable()->GetNucleusMass(theZ,theA);
    … …  
    107104                                      (U-delta0));
    108105                                                                 
    109     // compute the integrated probability of evaporation channel
     106
    110107    G4double RN = 1.5*fermi;
    111108
    … …  
    133130       
    134131    return Width;
    135 }
    136 
    137 
    138 
    139 
     132             
     133 } else if (OPTxs==1 || OPTxs==2 ||OPTxs==3 || OPTxs==4) {
     134
     135   G4double limit;
     136   if (useSICB) limit=theCoulombBarrierptr->GetCoulombBarrier(G4lrint(ResidualA),G4lrint(ResidualZ),U);
     137   else limit=0.;
     138
     139   if (MaximalKineticEnergy <= limit)  return 0.0;
     140
     141
     142   // if Coulomb barrier cutoff is superimposed for all cross sections the limit is the Coulomb Barrier
     143   G4double LowerLimit= limit;
     144
     145   //  MaximalKineticEnergy: asimptotic value (already accounted for in G4EvaporationChannel)     
     146
     147   G4double UpperLimit = MaximalKineticEnergy;
     148
     149
     150   G4double Width = IntegrateEmissionProbability(fragment,LowerLimit,UpperLimit);
     151
     152   return Width;
     153 } else{
     154   std::ostringstream errOs;
     155   errOs << "Bad option for cross sections at evaporation"  <<G4endl;
     156   throw G4HadronicException(__FILE__, __LINE__, errOs.str());
     157 }
     158 
     159}
     160
     161/////////////////////////////////////////////////////////////////////
     162
     163G4double G4EvaporationProbability::
     164IntegrateEmissionProbability(const G4Fragment & fragment, const G4double & Low, const G4double & Up )
     165{
     166
     167  static const G4int N = 10;
     168  // 10-Points Gauss-Legendre abcisas and weights
     169  static const G4double w[N] = {
     170    0.0666713443086881,
     171    0.149451349150581,
     172    0.219086362515982,
     173    0.269266719309996,
     174    0.295524224714753,
     175    0.295524224714753,
     176    0.269266719309996,
     177    0.219086362515982,
     178    0.149451349150581,
     179    0.0666713443086881
     180  };
     181  static const G4double x[N] = {
     182    -0.973906528517172,
     183    -0.865063366688985,
     184    -0.679409568299024,
     185    -0.433395394129247,
     186    -0.148874338981631,
     187    0.148874338981631,
     188    0.433395394129247,
     189    0.679409568299024,
     190    0.865063366688985,
     191    0.973906528517172
     192  };
     193
     194  G4double Total = 0.0;
     195
     196
     197  for (G4int i = 0; i < N; i++)
     198    {
     199
     200      G4double KineticE = ((Up-Low)*x[i]+(Up+Low))/2.0;
     201
     202      Total += w[i]*ProbabilityDistributionFunction(fragment, KineticE);
     203
     204    }
     205  Total *= (Up-Low)/2.0;
     206  return Total;
     207}
     208
     209
     210/////////////////////////////////////////////////////////
     211//New method (OPT=1,2,3,4)
     212
     213G4double G4EvaporationProbability::ProbabilityDistributionFunction( const G4Fragment & fragment, const G4double K)
     214{
     215
     216 
     217
     218
     219  G4double ResidualA = fragment.GetA() - theA;
     220  G4double ResidualZ = fragment.GetZ() - theZ; 
     221  G4double U = fragment.GetExcitationEnergy();
     222
     223  //        if(K <= theCoulombBarrierptr->GetCoulombBarrier(G4lrint(ResidualA),G4lrint(ResidualZ),U)) return 0.0;   
     224
     225  G4double delta1 = G4PairingCorrection::GetInstance()->GetPairingCorrection(static_cast<G4int>(ResidualA),static_cast<G4int>(ResidualZ));
     226
     227 
     228  G4double delta0 = G4PairingCorrection::GetInstance()->GetPairingCorrection(static_cast<G4int>(fragment.GetA()),static_cast<G4int>(fragment.GetZ()));
     229
     230 
     231  G4double ParticleMass = G4ParticleTable::GetParticleTable()->GetIonTable()->GetNucleusMass(theZ,theA);
     232
     233  G4double theSeparationEnergy= G4NucleiProperties::GetMassExcess(static_cast<G4int>(theA),static_cast<G4int>(theZ)) +
     234    G4NucleiProperties::GetMassExcess(static_cast<G4int>(ResidualA),static_cast<G4int>(ResidualZ)) -
     235    G4NucleiProperties::GetMassExcess(static_cast<G4int>(fragment.GetA()),static_cast<G4int>(fragment.GetZ()));
     236
     237  G4double a0 = theEvapLDPptr->LevelDensityParameter(static_cast<G4int>(fragment.GetA()),static_cast<G4int>(fragment.GetZ()),U - delta0);
     238
     239  G4double a1 = theEvapLDPptr->LevelDensityParameter(static_cast<G4int>(ResidualA),static_cast<G4int>(ResidualZ),U - theSeparationEnergy - delta1);
     240 
     241 
     242  G4double E0=U-delta0;
     243
     244  G4double E1=U-theSeparationEnergy-delta1-K;
     245
     246  if (E1<0.) return 0.;
     247
     248  //JMQ 14/02/09 BUG fixed: hbarc should be in the denominator instead of hbar_Planck
     249  //Without 1/hbar_Panck remains as a width
     250  //  G4double  Prob=Gamma*ParticleMass/((pi*hbar_Planck)*(pi*hbar_Planck)*
     251  //std::exp(2*std::sqrt(a0*E0)))*K*CrossSection(fragment,K)*std::exp(2*std::sqrt(a1*E1))*millibarn;
     252
     253  G4double Prob=Gamma*ParticleMass/((pi*hbarc)*(pi*hbarc)*std::exp(2*std::sqrt(a0*E0)))
     254    *K*CrossSection(fragment,K)*std::exp(2*std::sqrt(a1*E1))*millibarn;
     255
     256  return Prob;
     257}
     258
     259

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

    r819 r962  
    2626//
    2727// $Id: G4He3EvaporationChannel.cc,v 1.4 2006/06/29 20:10:33 gunter Exp $
    28 // GEANT4 tag $Name: $
     28// GEANT4 tag $Name: geant4-09-02-ref-02 $
    2929//
    3030// Hadronic Process: Nuclear De-excitations

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

    r819 r962  
    2424// ********************************************************************
    2525//
    26 //
    27 // $Id: G4He3EvaporationProbability.cc,v 1.4 2006/06/29 20:10:35 gunter Exp $
    28 // GEANT4 tag $Name: geant4-09-01-patch-02 $
     26//J.M. Quesada (August2008). Based on:
    2927//
    3028// Hadronic Process: Nuclear De-excitations
    31 // by V. Lara (Nov 1999)
    32 //
    33 
     29// by V. Lara (Oct 1998)
     30//
     31// Modif (03 September 2008) by J. M. Quesada for external choice of inverse
     32// cross section option
    3433
    3534#include "G4He3EvaporationProbability.hh"
    3635
    3736G4He3EvaporationProbability::G4He3EvaporationProbability() :
    38     G4EvaporationProbability(3,2,6) // A,Z,Gamma
    39 {
    40     std::vector<G4double>::size_type NumExcitedStatesEnergy = 31+1;
    41     std::vector<G4int>::size_type NumExcitedStatesSpin = 31+1;
    42     ExcitEnergies.reserve(NumExcitedStatesEnergy);
    43     ExcitSpins.reserve(NumExcitedStatesSpin);
    44     ExcitEnergies.insert(ExcitEnergies.begin(),NumExcitedStatesEnergy,0.0);
    45     ExcitSpins.insert(ExcitSpins.begin(),NumExcitedStatesSpin,0);
    46 
    47 
    48        
    49     ExcitEnergies[18] = 7.29*MeV;
    50     ExcitEnergies[20] = 6.48*MeV;
    51     ExcitEnergies[25] = 5.69*MeV;
    52     ExcitEnergies[26] = 8.31*MeV;
    53     ExcitEnergies[31] = 5.10*MeV;
    54 
    55     ExcitSpins[18] = 6;
    56     ExcitSpins[20] = 8;
    57     ExcitSpins[25] = 3;
    58     ExcitSpins[26] = 2;
    59     ExcitSpins[31] = 7;
    60        
    61     SetExcitationEnergiesPtr(&ExcitEnergies);
    62     SetExcitationSpinsPtr(&ExcitSpins);
     37   G4EvaporationProbability(3,2,2,&theCoulombBarrier) // A,Z,Gamma,&theCoulombBarrier
     38{
     39
    6340}
    6441
    … …  
    8966}
    9067
    91 G4double G4He3EvaporationProbability::CCoeficient(const G4double aZ) const
     68  G4double G4He3EvaporationProbability::CalcAlphaParam(const G4Fragment & fragment) 
     69  { return 1.0 + CCoeficient(static_cast<G4double>(fragment.GetZ()-GetZ()));}
     70       
     71  G4double G4He3EvaporationProbability::CalcBetaParam(const G4Fragment & ) 
     72  { return 0.0; }
     73
     74
     75G4double G4He3EvaporationProbability::CCoeficient(const G4double aZ)
    9276{
    9377    // Data comes from
    … …  
    11397    return C*(4.0/3.0);
    11498}
     99
     100///////////////////////////////////////////////////////////////////////////////////
     101//J. M. Quesada (Dec 2007-June 2008): New inverse reaction cross sections
     102//OPT=0 Dostrovski's parameterization
     103//OPT=1,2 Chatterjee's paramaterization
     104//OPT=3,4 Kalbach's parameterization
     105//
     106G4double 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 ************************
     134//(fitting to cross section from Bechetti & Greenles OM potential)
     135
     136G4double G4He3EvaporationProbability::GetOpt12(const  G4double K)
     137{
     138
     139G4double 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;}
     171             
     172   return xs;
     173
     174}
     175
     176// *********** OPT=3,4 : Kalbach's cross sections (from PRECO code)*************
     177G4double G4He3EvaporationProbability::GetOpt34(const  G4double K)
     178//c     ** 3he from o.m. of gibson et al
     179{
     180 
     181  G4double landa, mu, nu, p , signor(1.),sig;
     182  G4double ec,ecsq,xnulam,etest(0.),a;
     183  G4double b,ecut,cut,ecut2,geom,elab;
     184 
     185 
     186  G4double     flow = 1.e-18;
     187  G4double     spill= 1.e+18;
     188
     189
     190  G4double     p0 = -2.88;
     191  G4double     p1 = 205.6;
     192  G4double     p2 = -1487.;
     193  G4double     landa0 = 0.00459;
     194  G4double     landa1 = -8.93;
     195  G4double     mu0 = 611.2;
     196  G4double     mu1 = 0.35;
     197  G4double     nu0 = 473.8;
     198  G4double     nu1 = -468.2;
     199  G4double     nu2 = -2.225;     
     200 
     201  G4double      ra=0.80;
     202 
     203  //JMQ 13/02/09 increase of reduced radius to lower the barrier
     204  // ec = 1.44 * theZ * ResidualZ / (1.5*ResidualAthrd+ra);
     205  ec = 1.44 * theZ * ResidualZ / (1.7*ResidualAthrd+ra);
     206  ecsq = ec * ec;
     207  p = p0 + p1/ec + p2/ecsq;
     208  landa = landa0*ResidualA + landa1;
     209  a = std::pow(ResidualA,mu1);
     210  mu = mu0 * a;
     211  nu = a* (nu0+nu1*ec+nu2*ecsq); 
     212  xnulam = nu / landa;
     213  if (xnulam > spill) xnulam=0.;
     214  if (xnulam >= flow) etest = 1.2 *std::sqrt(xnulam);
     215 
     216  a = -2.*p*ec + landa - nu/ecsq;
     217  b = p*ecsq + mu + 2.*nu/ec;
     218  ecut = 0.;
     219  cut = a*a - 4.*p*b;
     220  if (cut > 0.) ecut = std::sqrt(cut);
     221  ecut = (ecut-a) / (p+p);
     222  ecut2 = ecut;
     223  if (cut < 0.) ecut2 = ecut - 2.;
     224  elab = K * FragmentA / ResidualA;
     225  sig = 0.;
     226
     227  if (elab <= ec) { //start for E<Ec
     228    if (elab > ecut2)  sig = (p*elab*elab+a*elab+b) * signor;
     229  }           //end for E<Ec
     230  else {           //start for E>Ec
     231    sig = (landa*elab+mu+nu/elab) * signor;
     232    geom = 0.;
     233    if (xnulam < flow || elab < etest) return sig;
     234    geom = std::sqrt(theA*K);
     235    geom = 1.23*ResidualAthrd + ra + 4.573/geom;
     236    geom = 31.416 * geom * geom;
     237    sig = std::max(geom,sig);
     238  }           //end for E>Ec
     239  return sig;
     240 
     241}
     242
     243//   ************************** end of cross sections *******************************
     244

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

    r819 r962  
    2626//
    2727// $Id: G4NeutronEvaporationChannel.cc,v 1.4 2006/06/29 20:10:37 gunter Exp $
    28 // GEANT4 tag $Name: $
     28// GEANT4 tag $Name: geant4-09-02-ref-02 $
    2929//
    3030// Hadronic Process: Nuclear De-excitations

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

    r819 r962  
    2424// ********************************************************************
    2525//
    26 //
    27 // $Id: G4NeutronEvaporationProbability.cc,v 1.4 2006/06/29 20:10:39 gunter Exp $
    28 // GEANT4 tag $Name: geant4-09-01-patch-02 $
     26//J.M. Quesada (August2008). Based on:
    2927//
    3028// Hadronic Process: Nuclear De-excitations
    31 // by V. Lara (Nov 1999)
    32 //
    33 
     29// by V. Lara (Oct 1998)
     30//
     31// Modif (03 September 2008) by J. M. Quesada for external choice of inverse
     32// cross section option
    3433
    3534#include "G4NeutronEvaporationProbability.hh"
    3635
    3736G4NeutronEvaporationProbability::G4NeutronEvaporationProbability() :
    38     G4EvaporationProbability(1,0,2) // A,Z,Gamma
    39 {
    40     std::vector<G4double>::size_type NumExcitedStatesEnergy = 31+1;
    41     std::vector<G4int>::size_type NumExcitedStatesSpin = 31+1;
    42     ExcitEnergies.reserve(NumExcitedStatesEnergy);
    43     ExcitSpins.reserve(NumExcitedStatesSpin);
    44     ExcitEnergies.insert(ExcitEnergies.begin(),NumExcitedStatesEnergy,0.0);
    45     ExcitSpins.insert(ExcitSpins.begin(),NumExcitedStatesSpin,0);
    46 
    47        
    48     ExcitEnergies[ 9] = 3.56*MeV;
    49     ExcitEnergies[10] = 0.48*MeV;
    50     ExcitEnergies[11] = 0.98*MeV;
    51     ExcitEnergies[12] = 0.43*MeV;
    52     ExcitEnergies[15] = 3.37*MeV;
    53     ExcitEnergies[17] = 0.72*MeV;
    54     ExcitEnergies[18] = 2.13*MeV;
    55     ExcitEnergies[19] = 0.95*MeV;
    56     ExcitEnergies[20] = 2.00*MeV;
    57     ExcitEnergies[21] = 4.44*MeV;
    58     ExcitEnergies[22] = 3.09*MeV;
    59     ExcitEnergies[23] = 6.09*MeV;
    60     ExcitEnergies[25] = 2.31*MeV;
    61     ExcitEnergies[26] = 5.28*MeV;
    62     ExcitEnergies[27] = 0.12*MeV;
    63     ExcitEnergies[28] = 5.22*MeV;
    64     ExcitEnergies[29] = 6.10*MeV;
    65     ExcitEnergies[30] = 0.87*MeV;
    66     ExcitEnergies[31] = 1.98*MeV;
    67 
    68 
    69     ExcitSpins[ 9] = 1;
    70     ExcitSpins[10] = 2;
    71     ExcitSpins[11] = 3;
    72     ExcitSpins[12] = 2;
    73     ExcitSpins[15] = 5;
    74     ExcitSpins[17] = 3;
    75     ExcitSpins[18] = 2;
    76     ExcitSpins[19] = 5;
    77     ExcitSpins[20] = 2;
    78     ExcitSpins[21] = 5;
    79     ExcitSpins[22] = 2;
    80     ExcitSpins[23] = 3;
    81     ExcitSpins[25] = 1;
    82     ExcitSpins[26] = 8;
    83     ExcitSpins[27] = 1;
    84     ExcitSpins[28] = 8;
    85     ExcitSpins[29] = 8;
    86     ExcitSpins[30] = 2;
    87     ExcitSpins[31] = 5;
    88 
    89     SetExcitationEnergiesPtr(&ExcitEnergies);
    90     SetExcitationSpinsPtr(&ExcitSpins);
     37    G4EvaporationProbability(1,0,2,&theCoulombBarrier) // A,Z,Gamma,&theCoulombBarrier
     38{
     39 
    9140}
    9241
    … …  
    11766    return true;
    11867}
     68
     69 G4double G4NeutronEvaporationProbability::CalcAlphaParam(const G4Fragment & fragment)
     70  { return 0.76+2.2/std::pow(static_cast<G4double>(fragment.GetA()-GetA()),1.0/3.0);}
     71
     72       
     73  G4double G4NeutronEvaporationProbability::CalcBetaParam(const G4Fragment &  fragment)
     74  { return (2.12/std::pow(static_cast<G4double>(fragment.GetA()-GetA()),2.0/3.0) - 0.05)*MeV/
     75      CalcAlphaParam(fragment); }
     76
     77
     78////////////////////////////////////////////////////////////////////////////////////
     79//J. M. Quesada (Dec 2007-June 2008): New inverse reaction cross sections
     80//OPT=0 Dostrovski's parameterization
     81//OPT=1,2 Chatterjee's paramaterization
     82//OPT=3,4 Kalbach's parameterization
     83//
     84 G4double G4NeutronEvaporationProbability::CrossSection(const  G4Fragment & fragment, const  G4double K)
     85{
     86  theA=GetA();
     87  theZ=GetZ();
     88  ResidualA=fragment.GetA()-theA;
     89  ResidualZ=fragment.GetZ()-theZ;
     90 
     91  ResidualAthrd=std::pow(ResidualA,0.33333);
     92  FragmentA=fragment.GetA();
     93  FragmentAthrd=std::pow(FragmentA,0.33333);
     94
     95
     96  if (OPTxs==0) {std::ostringstream errOs;
     97    errOs << "We should'n be here (OPT =0) at evaporation cross section calculation (neutrons)!!"  <<G4endl;
     98    throw G4HadronicException(__FILE__, __LINE__, errOs.str());
     99    return 0.;}
     100  else if( OPTxs==1 ||OPTxs==2) return GetOpt12( K);
     101  else if (OPTxs==3 || OPTxs==4)  return GetOpt34( K);
     102  else{
     103    std::ostringstream errOs;
     104    errOs << "BAD NEUTRON CROSS SECTION OPTION AT EVAPORATION!!"  <<G4endl;
     105    throw G4HadronicException(__FILE__, __LINE__, errOs.str());
     106    return 0.;
     107  }
     108}
     109
     110
     111
     112
     113
     114
     115//********************* OPT=1,2 : Chatterjee's cross section ************************
     116//(fitting to cross section from Bechetti & Greenles OM potential)
     117
     118G4double G4NeutronEvaporationProbability::GetOpt12(const  G4double K)
     119{
     120
     121  G4double Kc=K;
     122
     123// JMQ  xsec is set constat above limit of validity
     124  if (K>50) Kc=50;
     125
     126  G4double landa, landa0, landa1, mu, mu0, mu1,nu, nu0, nu1, nu2,xs;
     127
     128  landa0 = 18.57;
     129  landa1 = -22.93;
     130  mu0 = 381.7;
     131  mu1 = 24.31;
     132  nu0 = 0.172;
     133  nu1 = -15.39;
     134  nu2 = 804.8;
     135  landa = landa0/ResidualAthrd + landa1;
     136  mu = mu0*ResidualAthrd + mu1*ResidualAthrd*ResidualAthrd;
     137  nu = nu0*ResidualAthrd*ResidualA + nu1*ResidualAthrd*ResidualAthrd + nu2 ;
     138  xs=landa*Kc + mu + nu/Kc;
     139  if (xs <= 0.0 ){
     140    std::ostringstream errOs;
     141    G4cout<<"WARNING:  NEGATIVE OPT=1 neutron cross section "<<G4endl;     
     142    errOs << "RESIDUAL: Ar=" << ResidualA << " Zr=" << ResidualZ <<G4endl;
     143    errOs <<"  xsec("<<Kc<<" MeV) ="<<xs <<G4endl;
     144    throw G4HadronicException(__FILE__, __LINE__, errOs.str());
     145  }
     146  return xs;
     147}
     148
     149
     150// *********** OPT=3,4 : Kalbach's cross sections (from PRECO code)*************
     151G4double G4NeutronEvaporationProbability::GetOpt34(const  G4double K)
     152{
     153
     154  G4double landa, landa0, landa1, mu, mu0, mu1,nu, nu0, nu1, nu2;
     155  G4double p, p0, p1, p2;
     156  G4double flow,spill,ec,ecsq,xnulam,etest(0.),ra(0.),a,signor(1.),sig;
     157  G4double b,ecut,cut,ecut2,geom,elab;
     158
     159//safety initialization
     160  landa0=0;
     161  landa1=0;
     162  mu0=0.;
     163  mu1=0.;
     164  nu0=0.;
     165  nu1=0.;
     166  nu2=0.;
     167  p0=0.;
     168  p1=0.;
     169  p2=0.;
     170
     171
     172  flow = 1.e-18;
     173  spill= 1.0e+18;
     174
     175 
     176
     177// PRECO xs for neutrons is choosen
     178
     179  p0 = -312.;
     180  p1= 0.;
     181  p2 = 0.;
     182  landa0 = 12.10;
     183  landa1=  -11.27;
     184  mu0 = 234.1;
     185  mu1 = 38.26;
     186  nu0 = 1.55;
     187  nu1 = -106.1;
     188  nu2 = 1280.8;
     189
     190  if (ResidualA < 40.) signor=0.7+ResidualA*0.0075;
     191  if (ResidualA > 210.) signor = 1. + (ResidualA-210.)/250.;
     192  landa = landa0/ResidualAthrd + landa1;
     193  mu = mu0*ResidualAthrd + mu1*ResidualAthrd*ResidualAthrd;
     194  nu = nu0*ResidualAthrd*ResidualA + nu1*ResidualAthrd*ResidualAthrd + nu2;
     195
     196  // JMQ very low energy behaviour corrected (problem  for A (apprx.)>60)
     197  if (nu < 0.)nu=-nu;
     198
     199  ec = 0.5;
     200  ecsq = 0.25;
     201  p = p0;
     202  xnulam = 1.;
     203  etest = 32.;
     204//          ** etest is the energy above which the rxn cross section is
     205//          ** compared with the geometrical limit and the max taken.
     206//          ** xnulam here is a dummy value to be used later.
     207
     208
     209
     210  a = -2.*p*ec + landa - nu/ecsq;
     211  b = p*ecsq + mu + 2.*nu/ec;
     212  ecut = 0.;
     213  cut = a*a - 4.*p*b;
     214  if (cut > 0.) ecut = std::sqrt(cut);
     215  ecut = (ecut-a) / (p+p);
     216  ecut2 = ecut;
     217  if (cut < 0.) ecut2 = ecut - 2.;
     218  elab = K * FragmentA / ResidualA;
     219  sig = 0.;
     220  if (elab <= ec) { //start for E<Ec
     221    if (elab > ecut2)  sig = (p*elab*elab+a*elab+b) * signor;   
     222  }              //end for E<Ec
     223  else {           //start for  E>Ec
     224    sig = (landa*elab+mu+nu/elab) * signor;
     225    geom = 0.;
     226    if (xnulam < flow || elab < etest) return sig;
     227    geom = std::sqrt(theA*K);
     228    geom = 1.23*ResidualAthrd + ra + 4.573/geom;
     229    geom = 31.416 * geom * geom;
     230    sig = std::max(geom,sig);
     231  }
     232  return sig;}
     233
     234//   ************************** end of cross sections *******************************
     235

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

    r819 r962  
    2626//
    2727// $Id: G4ProtonEvaporationChannel.cc,v 1.4 2006/06/29 20:10:41 gunter Exp $
    28 // GEANT4 tag $Name: $
     28// GEANT4 tag $Name: geant4-09-02-ref-02 $
    2929//
    3030// Hadronic Process: Nuclear De-excitations

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

    r819 r962  
    2424// ********************************************************************
    2525//
    26 //
    27 // $Id: G4ProtonEvaporationProbability.cc,v 1.4 2006/06/29 20:10:43 gunter Exp $
    28 // GEANT4 tag $Name: geant4-09-01-patch-02 $
     26//J.M. Quesada (August2008). Based on:
    2927//
    3028// Hadronic Process: Nuclear De-excitations
    31 // by V. Lara (Nov 1999)
    32 //
    33 
     29// by V. Lara (Oct 1998)
     30//
     31// Modif (03 September 2008) by J. M. Quesada for external choice of inverse
     32// cross section option
    3433
    3534#include "G4ProtonEvaporationProbability.hh"
    3635
    3736G4ProtonEvaporationProbability::G4ProtonEvaporationProbability() :
    38     G4EvaporationProbability(1,1,2) // A,Z,Gamma
    39 {
    40     std::vector<G4double>::size_type NumExcitedStatesEnergy = 31+1;
    41     std::vector<G4int>::size_type NumExcitedStatesSpin = 31+1;
    42     ExcitEnergies.reserve(NumExcitedStatesEnergy);
    43     ExcitSpins.reserve(NumExcitedStatesSpin);
    44     ExcitEnergies.insert(ExcitEnergies.begin(),NumExcitedStatesEnergy,0.0);
    45     ExcitSpins.insert(ExcitSpins.begin(),NumExcitedStatesSpin,0);
    46 
    47 
    48 
    49     ExcitEnergies[15] = 5.96*MeV;
    50     ExcitEnergies[17] = 1.74*MeV;
    51     ExcitEnergies[18] = 4.44*MeV;
    52     ExcitEnergies[19] = 1.67*MeV;
    53     ExcitEnergies[20] = 4.32*MeV;
    54     ExcitEnergies[22] = 3.68*MeV;
    55     ExcitEnergies[23] = 6.69*MeV;
    56     ExcitEnergies[25] = 3.95*MeV;
    57     ExcitEnergies[26] = 6.32*MeV;
    58     ExcitEnergies[27] = 0.30*MeV;
    59     ExcitEnergies[28] = 6.18*MeV;
    60     ExcitEnergies[29] = 6.92*MeV;
    61     ExcitEnergies[30] = 3.06*MeV;
    62     ExcitEnergies[31] = 3.57*MeV;
    63 
    64 
    65     ExcitSpins[15] = 8;
    66     ExcitSpins[17] = 1;
    67     ExcitSpins[18] = 6;
    68     ExcitSpins[19] = 5;
    69     ExcitSpins[20] = 6;
    70     ExcitSpins[22] = 4;
    71     ExcitSpins[23] = 8;
    72     ExcitSpins[25] = 3;
    73     ExcitSpins[26] = 4;
    74     ExcitSpins[27] = 7;
    75     ExcitSpins[28] = 4;
    76     ExcitSpins[29] = 5;
    77     ExcitSpins[30] = 2;
    78     ExcitSpins[31] = 10;
    79        
    80     SetExcitationEnergiesPtr(&ExcitEnergies);
    81     SetExcitationSpinsPtr(&ExcitSpins);
    82        
     37    G4EvaporationProbability(1,1,2,&theCoulombBarrier) // A,Z,Gamma,&theCoulombBarrier
     38{
     39   
    8340}
    8441
    … …  
    10663}
    10764
    108 G4double G4ProtonEvaporationProbability::CCoeficient(const G4double aZ) const
     65  G4double G4ProtonEvaporationProbability::CalcAlphaParam(const G4Fragment & fragment)
     66  { return 1.0 + CCoeficient(static_cast<G4double>(fragment.GetZ()-GetZ()));}
     67       
     68  G4double G4ProtonEvaporationProbability::CalcBetaParam(const G4Fragment & ) 
     69  { return 0.0; }
     70
     71  G4double G4ProtonEvaporationProbability::CCoeficient(const G4double aZ)
    10972{
    11073    // Data comes from
    … …  
    12689       
    12790}
     91
     92///////////////////////////////////////////////////////////////////////////////////
     93//J. M. Quesada (Dec 2007-June 2008): New inverse reaction cross sections
     94//OPT=0 Dostrovski's parameterization
     95//OPT=1,2 Chatterjee's paramaterization
     96//OPT=3,4 Kalbach's parameterization
     97//
     98G4double G4ProtonEvaporationProbability::CrossSection(const  G4Fragment & fragment, const  G4double K)
     99{
     100//  G4cout<<" In G4ProtonEVaporationProbability OPTxs="<<OPTxs<<G4endl;
     101//  G4cout<<" In G4ProtonEVaporationProbability useSICB="<<useSICB<<G4endl;
     102
     103  theA=GetA();
     104  theZ=GetZ();
     105  ResidualA=fragment.GetA()-theA;
     106  ResidualZ=fragment.GetZ()-theZ;
     107 
     108  ResidualAthrd=std::pow(ResidualA,0.33333);
     109  FragmentA=fragment.GetA();
     110  FragmentAthrd=std::pow(FragmentA,0.33333);
     111  U=fragment.GetExcitationEnergy();
     112
     113
     114  if (OPTxs==0) {std::ostringstream errOs;
     115    errOs << "We should'n be here (OPT =0) at evaporation cross section calculation (protons)!!"  <<G4endl;
     116    throw G4HadronicException(__FILE__, __LINE__, errOs.str());
     117    return 0.;}
     118  else if( OPTxs==1 ) return GetOpt1( K);
     119  else if( OPTxs==2 ||OPTxs==4) return GetOpt2( K);
     120  else if (OPTxs==3 )  return GetOpt3( K);
     121  else{
     122    std::ostringstream errOs;
     123    errOs << "BAD PROTON CROSS SECTION OPTION AT EVAPORATION!!"  <<G4endl;
     124    throw G4HadronicException(__FILE__, __LINE__, errOs.str());
     125    return 0.;
     126  }
     127}
     128//********************* OPT=1 : Chatterjee's cross section ************************
     129//(fitting to cross section from Bechetti & Greenles OM potential)
     130
     131G4double G4ProtonEvaporationProbability::GetOpt1(const  G4double K)
     132{
     133  G4double Kc=K;
     134
     135// JMQ  xsec is set constat above limit of validity
     136  if (K>50)  Kc=50;
     137
     138  G4double landa, landa0, landa1, mu, mu0, mu1,nu, nu0, nu1, nu2,xs;
     139  G4double p, p0, p1, p2,Ec,delta,q,r,ji;
     140 
     141  p0 = 15.72;
     142  p1 = 9.65;
     143  p2 = -449.0;
     144  landa0 = 0.00437;
     145  landa1 = -16.58;
     146  mu0 = 244.7;
     147  mu1 = 0.503;
     148  nu0 = 273.1;
     149  nu1 = -182.4;
     150  nu2 = -1.872; 
     151  delta=0.; 
     152
     153  Ec = 1.44*theZ*ResidualZ/(1.5*ResidualAthrd+delta);
     154  p = p0 + p1/Ec + p2/(Ec*Ec);
     155  landa = landa0*ResidualA + landa1;
     156  mu = mu0*std::pow(ResidualA,mu1);
     157  nu = std::pow(ResidualA,mu1)*(nu0 + nu1*Ec + nu2*(Ec*Ec));
     158  q = landa - nu/(Ec*Ec) - 2*p*Ec;
     159  r = mu + 2*nu/Ec + p*(Ec*Ec);
     160
     161  ji=std::max(Kc,Ec);
     162  if(Kc < Ec) { xs = p*Kc*Kc + q*Kc + r;}
     163  else {xs = p*(Kc - ji)*(Kc - ji) + landa*Kc + mu + nu*(2 - Kc/ji)/ji ;}
     164   if (xs <0.0) {xs=0.0;}
     165
     166   return xs;
     167
     168}
     169
     170
     171
     172//************* OPT=2 : Wellisch's proton reaction cross section ************************
     173
     174G4double G4ProtonEvaporationProbability::GetOpt2(const  G4double K)
     175{
     176
     177  G4double rnpro,rnneu,eekin,ekin,ff1,ff2,ff3,r0,fac,fac1,fac2,b0,xine_th(0);
     178 
     179//This is redundant when the Coulomb  barrier is overimposed to all cross sections
     180//It should be kept when Coulomb barrier only imposed at OPTxs=2, this is why ..
     181
     182         if(!useSICB && K <= theCoulombBarrier.GetCoulombBarrier(G4lrint(ResidualA),G4lrint(ResidualZ),U)) return xine_th=0.0;
     183
     184  eekin=K;
     185  rnpro=ResidualZ;
     186  rnneu=ResidualA-ResidualZ;
     187  ekin=eekin/1000;
     188
     189  r0=1.36*1.e-15;
     190  fac=pi*r0*r0;
     191  b0=2.247-0.915*(1.-1./ResidualAthrd);
     192  fac1=b0*(1.-1./ResidualAthrd);
     193  fac2=1.;
     194  if(rnneu > 1.5) fac2=std::log(rnneu);
     195  xine_th= 1.e+31*fac*fac2*(1.+ResidualAthrd-fac1);
     196  xine_th=(1.-0.15*std::exp(-ekin))*xine_th/(1.00-0.0007*ResidualA);   
     197  ff1=0.70-0.0020*ResidualA ;
     198  ff2=1.00+1/ResidualA;
     199  ff3=0.8+18/ResidualA-0.002*ResidualA;
     200  fac=1.-(1./(1.+std::exp(-8.*ff1*(std::log10(ekin)+1.37*ff2))));
     201  xine_th=xine_th*(1.+ff3*fac);
     202  ff1=1.-1/ResidualA-0.001*ResidualA;
     203  ff2=1.17-2.7/ResidualA-0.0014*ResidualA;
     204  fac=-8.*ff1*(std::log10(ekin)+2.0*ff2);
     205  fac=1./(1.+std::exp(fac));
     206  xine_th=xine_th*fac;               
     207  if (xine_th < 0.0){
     208    std::ostringstream errOs;
     209    G4cout<<"WARNING:  negative Wellisch cross section "<<G4endl;
     210    errOs << "RESIDUAL: A=" << ResidualA << " Z=" << ResidualZ <<G4endl;
     211    errOs <<"  xsec("<<ekin<<" MeV) ="<<xine_th <<G4endl;
     212    throw G4HadronicException(__FILE__, __LINE__, errOs.str());
     213  }
     214
     215  return xine_th;
     216           
     217}
     218
     219
     220// *********** OPT=3 : Kalbach's cross sections (from PRECO code)*************
     221G4double G4ProtonEvaporationProbability::GetOpt3(const  G4double K)
     222{
     223//     ** p from  becchetti and greenlees (but modified with sub-barrier
     224//     ** correction function and xp2 changed from -449)
     225// JMQ (june 2008) : refinement of proton cross section for light systems
     226//
     227  G4double landa, landa0, landa1, mu, mu0, mu1,nu, nu0, nu1, nu2;
     228  G4double p, p0, p1, p2;
     229  p0 = 15.72;
     230  p1 = 9.65;
     231  p2 = -300.;
     232  landa0 = 0.00437;
     233  landa1 = -16.58;
     234  mu0 = 244.7;
     235  mu1 = 0.503;
     236  nu0 = 273.1;
     237  nu1 = -182.4;
     238  nu2 = -1.872;
     239
     240// parameters for  proton cross section refinement
     241  G4double afit,bfit,a2,b2;
     242  afit=-0.0785656;
     243  bfit=5.10789;
     244  a2= -0.00089076;
     245  b2= 0.0231597; 
     246
     247  G4double ec,ecsq,xnulam,etest(0.),ra(0.),a,w,c,signor(1.),signor2,sig;
     248  G4double b,ecut,cut,ecut2,geom,elab;
     249
     250
     251  G4double      flow = 1.e-18;
     252  G4double       spill= 1.e+18;
     253
     254
     255  if (ResidualA <= 60.)  signor = 0.92;
     256  else if (ResidualA < 100.) signor = 0.8 + ResidualA*0.002;
     257
     258
     259  ec = 1.44 * theZ * ResidualZ / (1.5*ResidualAthrd+ra);
     260  ecsq = ec * ec;
     261  p = p0 + p1/ec + p2/ecsq;
     262  landa = landa0*ResidualA + landa1;
     263  a = std::pow(ResidualA,mu1);
     264  mu = mu0 * a;
     265  nu = a* (nu0+nu1*ec+nu2*ecsq);
     266 
     267  c =std::min(3.15,ec*0.5);
     268  w = 0.7 * c / 3.15;
     269
     270  xnulam = nu / landa;
     271  if (xnulam > spill) xnulam=0.;
     272  if (xnulam >= flow) etest =std::sqrt(xnulam) + 7.;
     273
     274  a = -2.*p*ec + landa - nu/ecsq;
     275  b = p*ecsq + mu + 2.*nu/ec;
     276  ecut = 0.;
     277  cut = a*a - 4.*p*b;
     278  if (cut > 0.) ecut = std::sqrt(cut);
     279  ecut = (ecut-a) / (p+p);
     280  ecut2 = ecut;
     281  if (cut < 0.) ecut2 = ecut - 2.;
     282  elab = K * FragmentA / ResidualA;
     283  sig = 0.;
     284  if (elab <= ec) { //start for E<Ec
     285    if (elab > ecut2)  sig = (p*elab*elab+a*elab+b) * signor;
     286    signor2 = (ec-elab-c) / w;
     287    signor2 = 1. + std::exp(signor2);
     288    sig = sig / signor2;
     289    // signor2 is empirical global corr'n at low elab for protons in PRECO, not enough for light targets
     290    //  refinement for proton cross section
     291    if (ResidualZ<=26)
     292      sig = sig*std::exp(-(a2*ResidualZ + b2)*(elab-(afit*ResidualZ+bfit)*ec)*(elab-(afit*ResidualZ+bfit)*ec));     
     293                       }              //end for E<Ec
     294  else            {           //start for  E>Ec
     295    sig = (landa*elab+mu+nu/elab) * signor;
     296
     297    //  refinement for proton cross section
     298    if ( ResidualZ<=26 && elab <=(afit*ResidualZ+bfit)*ec)
     299           sig = sig*std::exp(-(a2*ResidualZ + b2)*(elab-(afit*ResidualZ+bfit)*ec)*(elab-(afit*ResidualZ+bfit)*ec));
     300                                                                               
     301//
     302    geom = 0.;
     303
     304    if (xnulam < flow || elab < etest)
     305     {
     306        if (sig <0.0) {sig=0.0;}
     307        return sig;
     308      }
     309    geom = std::sqrt(theA*K);
     310    geom = 1.23*ResidualAthrd + ra + 4.573/geom;
     311    geom = 31.416 * geom * geom;
     312    sig = std::max(geom,sig);
     313
     314  }   //end for E>Ec
     315 return sig;}
     316
     317
     318
     319//   ************************** end of cross sections *******************************
     320

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

    r819 r962  
    2626//
    2727// $Id: G4TritonEvaporationChannel.cc,v 1.4 2006/06/29 20:10:45 gunter Exp $
    28 // GEANT4 tag $Name: $
     28// GEANT4 tag $Name: geant4-09-02-ref-02 $
    2929//
    3030// Hadronic Process: Nuclear De-excitations

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

    r819 r962  
    2424// ********************************************************************
    2525//
    26 //
    27 // $Id: G4TritonEvaporationProbability.cc,v 1.4 2006/06/29 20:10:47 gunter Exp $
    28 // GEANT4 tag $Name: geant4-09-01-patch-02 $
     26//J.M. Quesada (August2008). Based on:
    2927//
    3028// Hadronic Process: Nuclear De-excitations
    31 // by V. Lara (Nov 1999)
    32 //
    33 
     29// by V. Lara (Oct 1998)
     30//
     31// Modif (03 September 2008) by J. M. Quesada for external choice of inverse
     32// cross section option
    3433
    3534#include "G4TritonEvaporationProbability.hh"
    3635
    3736G4TritonEvaporationProbability::G4TritonEvaporationProbability() :
    38     G4EvaporationProbability(3,1,6) // A,Z,Gamma
    39 {
    40     std::vector<G4double>::size_type NumExcitedStatesEnergy = 31+1;
    41     std::vector<G4int>::size_type NumExcitedStatesSpin = 31+1;
    42     ExcitEnergies.reserve(NumExcitedStatesEnergy);
    43     ExcitSpins.reserve(NumExcitedStatesSpin);
    44     ExcitEnergies.insert(ExcitEnergies.begin(),NumExcitedStatesEnergy,0.0);
    45     ExcitSpins.insert(ExcitSpins.begin(),NumExcitedStatesSpin,0);
    46 
    47        
    48     ExcitEnergies[15] = 6.26*MeV;
    49     ExcitEnergies[17] = 3.59*MeV;
    50     ExcitEnergies[18] = 6.76*MeV;
    51     ExcitEnergies[20] = 6.34*MeV;
    52     ExcitEnergies[23] = 7.34*MeV;
    53     ExcitEnergies[25] = 5.11*MeV;
    54     ExcitEnergies[26] = 7.57*MeV;
    55     ExcitEnergies[28] = 7.28*MeV;
    56     ExcitEnergies[31] = 4.46*MeV;
    57 
    58     ExcitSpins[15] = 5;
    59     ExcitSpins[17] = 5;
    60     ExcitSpins[18] = 10;
    61     ExcitSpins[20] = 2;
    62     ExcitSpins[23] = 5;
    63     ExcitSpins[25] = 5;
    64     ExcitSpins[26] = 8;
    65     ExcitSpins[28] = 8;
    66     ExcitSpins[31] = 3;
    67 
    68     SetExcitationEnergiesPtr(&ExcitEnergies);
    69     SetExcitationSpinsPtr(&ExcitSpins);
     37    G4EvaporationProbability(3,1,2,&theCoulombBarrier) // A,Z,Gamma,&theCoulombBarrier
     38{
     39
    7040}
    7141
    … …  
    9666}
    9767
    98 G4double G4TritonEvaporationProbability::CCoeficient(const G4double aZ) const
     68   G4double G4TritonEvaporationProbability::CalcAlphaParam(const G4Fragment & fragment)
     69  { return 1.0 + CCoeficient(static_cast<G4double>(fragment.GetZ()-GetZ()));}
     70       
     71  G4double G4TritonEvaporationProbability::CalcBetaParam(const G4Fragment & )
     72  { return 0.0; }
     73
     74G4double G4TritonEvaporationProbability::CCoeficient(const G4double aZ)
    9975{
    10076    // Data comes from
    … …  
    11793       
    11894}
     95
     96///////////////////////////////////////////////////////////////////////////////////
     97//J. M. Quesada (Dec 2007-June 2008): New inverse reaction cross sections
     98//OPT=0 Dostrovski's parameterization
     99//OPT=1,2 Chatterjee's paramaterization
     100//OPT=3,4 Kalbach's parameterization
     101//
     102G4double G4TritonEvaporationProbability::CrossSection(const  G4Fragment & fragment, const  G4double K)
     103{
     104  theA=GetA();
     105  theZ=GetZ();
     106  ResidualA=fragment.GetA()-theA;
     107  ResidualZ=fragment.GetZ()-theZ;
     108 
     109  ResidualAthrd=std::pow(ResidualA,0.33333);
     110  FragmentA=fragment.GetA();
     111  FragmentAthrd=std::pow(FragmentA,0.33333);
     112
     113 
     114  if (OPTxs==0) {std::ostringstream errOs;
     115    errOs << "We should'n be here (OPT =0) at evaporation cross section calculation (tritons)!!" 
     116          <<G4endl;
     117    throw G4HadronicException(__FILE__, __LINE__, errOs.str());
     118    return 0.;}
     119  if( OPTxs==1 || OPTxs==2) return G4TritonEvaporationProbability::GetOpt12( K);
     120  else if (OPTxs==3 || OPTxs==4)  return G4TritonEvaporationProbability::GetOpt34( K);
     121  else{
     122    std::ostringstream errOs;
     123    errOs << "BAD Triton CROSS SECTION OPTION AT EVAPORATION!!"  <<G4endl;
     124    throw G4HadronicException(__FILE__, __LINE__, errOs.str());
     125    return 0.;
     126  }
     127}
     128
     129//
     130//********************* OPT=1,2 : Chatterjee's cross section ************************
     131//(fitting to cross section from Bechetti & Greenles OM potential)
     132
     133G4double G4TritonEvaporationProbability::GetOpt12(const  G4double K)
     134{
     135
     136  G4double Kc=K;
     137
     138// JMQ xsec is set constat above limit of validity
     139  if (K>50) Kc=50;
     140
     141 G4double landa ,mu ,nu ,p , Ec,q,r,ji,xs;
     142//G4double Eo(0),epsilon1(0),epsilon2(0),discri(0);
     143
     144 
     145 G4double    p0 = -11.04;
     146 G4double    p1 = 619.1;
     147 G4double    p2 = -2147.;
     148 G4double    landa0 = -0.0426;
     149 G4double    landa1 = -10.33;
     150 G4double    mu0 = 601.9;
     151 G4double    mu1 = 0.37;
     152 G4double    nu0 = 583.0;
     153 G4double    nu1 = -546.2;
     154 G4double    nu2 = 1.718; 
     155 G4double    delta=1.2;           
     156
     157 Ec = 1.44*theZ*ResidualZ/(1.5*ResidualAthrd+delta);
     158 p = p0 + p1/Ec + p2/(Ec*Ec);
     159 landa = landa0*ResidualA + landa1;
     160 mu = mu0*std::pow(ResidualA,mu1);
     161 nu = std::pow(ResidualA,mu1)*(nu0 + nu1*Ec + nu2*(Ec*Ec));
     162 q = landa - nu/(Ec*Ec) - 2*p*Ec;
     163 r = mu + 2*nu/Ec + p*(Ec*Ec);
     164
     165 ji=std::max(Kc,Ec);
     166 if(Kc < Ec) { xs = p*Kc*Kc + q*Kc + r;}
     167 else {xs = p*(Kc - ji)*(Kc - ji) + landa*Kc + mu + nu*(2 - Kc/ji)/ji ;}
     168 
     169 if (xs <0.0) {xs=0.0;}
     170 
     171 return xs;
     172
     173}
     174
     175// *********** OPT=3,4 : Kalbach's cross sections (from PRECO code)*************
     176G4double G4TritonEvaporationProbability::GetOpt34(const  G4double K)
     177//     ** t from o.m. of hafele, flynn et al
     178{
     179 
     180  G4double landa, mu, nu, p , signor(1.),sig;
     181  G4double ec,ecsq,xnulam,etest(0.),a;
     182  G4double b,ecut,cut,ecut2,geom,elab;
     183 
     184 
     185  G4double     flow = 1.e-18;
     186  G4double     spill= 1.e+18;
     187
     188  G4double     p0 = -21.45;
     189  G4double     p1 = 484.7;
     190  G4double     p2 = -1608.;
     191  G4double     landa0 = 0.0186;
     192  G4double     landa1 = -8.90;
     193  G4double     mu0 = 686.3;
     194  G4double     mu1 = 0.325;
     195  G4double     nu0 = 368.9;
     196  G4double     nu1 = -522.2;
     197  G4double     nu2 = -4.998; 
     198 
     199  G4double      ra=0.80;
     200       
     201  //JMQ 13/02/09 increase of reduced radius to lower the barrier
     202  // ec = 1.44 * theZ * ResidualZ / (1.5*ResidualAthrd+ra);
     203  ec = 1.44 * theZ * ResidualZ / (1.7*ResidualAthrd+ra);
     204  ecsq = ec * ec;
     205  p = p0 + p1/ec + p2/ecsq;
     206  landa = landa0*ResidualA + landa1;
     207  a = std::pow(ResidualA,mu1);
     208  mu = mu0 * a;
     209  nu = a* (nu0+nu1*ec+nu2*ecsq); 
     210  xnulam = nu / landa;
     211  if (xnulam > spill) xnulam=0.;
     212  if (xnulam >= flow) etest = 1.2 *std::sqrt(xnulam);
     213 
     214  a = -2.*p*ec + landa - nu/ecsq;
     215  b = p*ecsq + mu + 2.*nu/ec;
     216  ecut = 0.;
     217  cut = a*a - 4.*p*b;
     218  if (cut > 0.) ecut = std::sqrt(cut);
     219  ecut = (ecut-a) / (p+p);
     220  ecut2 = ecut;
     221  if (cut < 0.) ecut2 = ecut - 2.;
     222  elab = K * FragmentA / ResidualA;
     223  sig = 0.;
     224
     225  if (elab <= ec) { //start for E<Ec
     226    if (elab > ecut2)  sig = (p*elab*elab+a*elab+b) * signor;
     227  }           //end for E<Ec
     228  else {           //start for E>Ec 
     229    sig = (landa*elab+mu+nu/elab) * signor;
     230    geom = 0.;
     231    if (xnulam < flow || elab < etest) return sig;
     232    geom = std::sqrt(theA*K);
     233    geom = 1.23*ResidualAthrd + ra + 4.573/geom;
     234    geom = 31.416 * geom * geom;
     235    sig = std::max(geom,sig);
     236  }           //end for E>Ec
     237  return sig;
     238 
     239}
     240
     241//   ************************** end of cross sections *******************************
     242

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

    r819 r962  
    2626//
    2727// $Id: G4VEvaporation.cc,v 1.5 2006/06/29 20:10:49 gunter Exp $
    28 // GEANT4 tag $Name: $
     28// GEANT4 tag $Name: geant4-09-02-ref-02 $
    2929//
    3030// Hadronic Process: Nuclear De-excitations
Note: See TracChangeset for help on using the changeset viewer.