Changeset 1347 for trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4NeutronEvaporationProbability.cc
- Timestamp:
- Dec 22, 2010, 3:52:27 PM (13 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4NeutronEvaporationProbability.cc
r962 r1347 24 24 // ******************************************************************** 25 25 // 26 //J.M. Quesada (August2008). Based on: 26 // $Id: G4NeutronEvaporationProbability.cc,v 1.16 2010/11/17 11:06:03 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-04-ref-00 $ 28 // 29 // J.M. Quesada (August2008). Based on: 27 30 // 28 31 // Hadronic Process: Nuclear De-excitations 29 32 // by V. Lara (Oct 1998) 30 33 // 31 // Modif (03 September 2008) by J. M. Quesada for external choice of inverse 32 // cross section option 34 // Modified: 35 // 03-09-2008 J.M. Quesada for external choice of inverse cross section option 36 // 17-11-2010 V.Ivanchenko integer Z and A 33 37 34 38 #include "G4NeutronEvaporationProbability.hh" … … 36 40 G4NeutronEvaporationProbability::G4NeutronEvaporationProbability() : 37 41 G4EvaporationProbability(1,0,2,&theCoulombBarrier) // A,Z,Gamma,&theCoulombBarrier 38 { 39 40 } 41 42 43 G4NeutronEvaporationProbability::G4NeutronEvaporationProbability(const G4NeutronEvaporationProbability &) : G4EvaporationProbability() 44 { 45 throw G4HadronicException(__FILE__, __LINE__, "G4NeutronEvaporationProbability::copy_constructor meant to not be accessable"); 46 } 47 48 49 50 51 const G4NeutronEvaporationProbability & G4NeutronEvaporationProbability:: 52 operator=(const G4NeutronEvaporationProbability &) 53 { 54 throw G4HadronicException(__FILE__, __LINE__, "G4NeutronEvaporationProbability::operator= meant to not be accessable"); 55 return *this; 56 } 57 58 59 G4bool G4NeutronEvaporationProbability::operator==(const G4NeutronEvaporationProbability &) const 60 { 61 return false; 62 } 63 64 G4bool G4NeutronEvaporationProbability::operator!=(const G4NeutronEvaporationProbability &) const 65 { 66 return true; 67 } 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 42 {} 43 44 G4NeutronEvaporationProbability::~G4NeutronEvaporationProbability() 45 {} 46 47 G4double G4NeutronEvaporationProbability::CalcAlphaParam(const G4Fragment & fragment) 48 { 49 return 0.76+2.2/fG4pow->Z13(fragment.GetA_asInt()-GetA()); 50 } 72 51 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 52 G4double G4NeutronEvaporationProbability::CalcBetaParam(const G4Fragment & fragment) 53 { 54 return (2.12/fG4pow->Z23(fragment.GetA_asInt()-GetA()) - 0.05)*MeV/ 55 CalcAlphaParam(fragment); 56 } 77 57 78 58 //////////////////////////////////////////////////////////////////////////////////// … … 82 62 //OPT=3,4 Kalbach's parameterization 83 63 // 84 G4double G4NeutronEvaporationProbability::CrossSection(const G4Fragment & fragment, const G4double K) 64 G4double 65 G4NeutronEvaporationProbability::CrossSection(const G4Fragment & fragment, G4double K) 85 66 { 86 67 theA=GetA(); 87 68 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 69 ResidualA=fragment.GetA_asInt()-theA; 70 ResidualZ=fragment.GetZ_asInt()-theZ; 71 72 ResidualAthrd=fG4pow->Z13(ResidualA); 73 FragmentA=fragment.GetA_asInt(); 74 FragmentAthrd=fG4pow->Z13(FragmentA); 95 75 96 76 if (OPTxs==0) {std::ostringstream errOs; … … 107 87 } 108 88 } 109 110 111 112 113 114 115 //********************* OPT=1,2 : Chatterjee's cross section ************************ 89 90 //********************* OPT=1,2 : Chatterjee's cross section *************** 116 91 //(fitting to cross section from Bechetti & Greenles OM potential) 117 92 118 G4double G4NeutronEvaporationProbability::GetOpt12( constG4double K)93 G4double G4NeutronEvaporationProbability::GetOpt12(G4double K) 119 94 { 120 121 95 G4double Kc=K; 122 96 123 // JMQ xsec is set constat above limit of validity 124 if (K>50) Kc=50; 97 // Pramana (Bechetti & Greenles) for neutrons is chosen 98 99 // JMQ xsec is set constat above limit of validity 100 if (K > 50*MeV) { Kc = 50*MeV; } 125 101 126 102 G4double landa, landa0, landa1, mu, mu0, mu1,nu, nu0, nu1, nu2,xs; … … 143 119 errOs <<" xsec("<<Kc<<" MeV) ="<<xs <<G4endl; 144 120 throw G4HadronicException(__FILE__, __LINE__, errOs.str()); 145 }121 } 146 122 return xs; 147 123 } 148 124 149 150 125 // *********** OPT=3,4 : Kalbach's cross sections (from PRECO code)************* 151 G4double G4NeutronEvaporationProbability::GetOpt34( constG4double K)126 G4double G4NeutronEvaporationProbability::GetOpt34(G4double K) 152 127 { 153 154 128 G4double landa, landa0, landa1, mu, mu0, mu1,nu, nu0, nu1, nu2; 155 129 G4double p, p0, p1, p2; … … 157 131 G4double b,ecut,cut,ecut2,geom,elab; 158 132 159 //safety initialization133 //safety initialization 160 134 landa0=0; 161 135 landa1=0; … … 169 143 p2=0.; 170 144 171 172 145 flow = 1.e-18; 173 146 spill= 1.0e+18; 174 147 175 176 177 // PRECO xs for neutrons is choosen 178 148 // PRECO xs for neutrons is choosen 179 149 p0 = -312.; 180 150 p1= 0.; … … 188 158 nu2 = 1280.8; 189 159 190 if (ResidualA < 40 .) signor=0.7+ResidualA*0.0075;191 if (ResidualA > 210 .) signor = 1. + (ResidualA-210.)/250.;160 if (ResidualA < 40) { signor =0.7 + ResidualA*0.0075; } 161 if (ResidualA > 210) { signor = 1. + (ResidualA-210)/250.; } 192 162 landa = landa0/ResidualAthrd + landa1; 193 163 mu = mu0*ResidualAthrd + mu1*ResidualAthrd*ResidualAthrd; … … 195 165 196 166 // JMQ very low energy behaviour corrected (problem for A (apprx.)>60) 197 if (nu < 0.) nu=-nu;167 if (nu < 0.) { nu=-nu; } 198 168 199 169 ec = 0.5; … … 202 172 xnulam = 1.; 203 173 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 174 // ** etest is the energy above which the rxn cross section is 175 // ** compared with the geometrical limit and the max taken. 176 // ** xnulam here is a dummy value to be used later. 209 177 210 178 a = -2.*p*ec + landa - nu/ecsq; … … 212 180 ecut = 0.; 213 181 cut = a*a - 4.*p*b; 214 if (cut > 0.) ecut = std::sqrt(cut);182 if (cut > 0.) { ecut = std::sqrt(cut); } 215 183 ecut = (ecut-a) / (p+p); 216 184 ecut2 = ecut; 217 if (cut < 0.) ecut2 = ecut - 2.;218 elab = K * FragmentA / ResidualA;185 if (cut < 0.) { ecut2 = ecut - 2.; } 186 elab = K * FragmentA / G4double(ResidualA); 219 187 sig = 0.; 220 188 if (elab <= ec) { //start for E<Ec 221 if (elab > ecut2) sig = (p*elab*elab+a*elab+b) * signor;189 if (elab > ecut2) { sig = (p*elab*elab+a*elab+b) * signor; } 222 190 } //end for E<Ec 223 191 else { //start for E>Ec 224 192 sig = (landa*elab+mu+nu/elab) * signor; 225 193 geom = 0.; 226 if (xnulam < flow || elab < etest) return sig;194 if (xnulam < flow || elab < etest) { return sig; } 227 195 geom = std::sqrt(theA*K); 228 196 geom = 1.23*ResidualAthrd + ra + 4.573/geom; 229 197 geom = 31.416 * geom * geom; 230 198 sig = std::max(geom,sig); 199 231 200 } 232 return sig;} 233 234 // ************************** end of cross sections ******************************* 235 201 return sig; 202 } 203
Note: See TracChangeset
for help on using the changeset viewer.