Changeset 1340 for trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/src/G4PreCompoundNeutron.cc
- Timestamp:
- Nov 5, 2010, 3:45:55 PM (14 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/src/G4PreCompoundNeutron.cc
r1337 r1340 24 24 // ******************************************************************** 25 25 // 26 // 27 // $Id: G4PreCompoundNeutron.cc,v 1.4 2009/02/11 18:06:00 vnivanch Exp $ 28 // GEANT4 tag $Name: geant4-09-04-beta-01 $ 26 // $Id: G4PreCompoundNeutron.cc,v 1.5 2010/08/28 15:16:55 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-03-ref-09 $ 29 28 // 30 29 // ------------------------------------------------------------------- … … 40 39 // 21.08.2008 J. M. Quesada add choice of options 41 40 // 10.02.2009 J. M. Quesada set default opt3 41 // 20.08.2010 V.Ivanchenko added G4Pow and G4PreCompoundParameters pointers 42 // use int Z and A and cleanup 42 43 // 43 44 44 45 #include "G4PreCompoundNeutron.hh" 45 46 G4ReactionProduct * G4PreCompoundNeutron::GetReactionProduct() const 47 { 48 G4ReactionProduct * theReactionProduct = 49 new G4ReactionProduct(G4Neutron::NeutronDefinition()); 50 theReactionProduct->SetMomentum(GetMomentum().vect()); 51 theReactionProduct->SetTotalEnergy(GetMomentum().e()); 52 #ifdef PRECOMPOUND_TEST 53 theReactionProduct->SetCreatorModel("G4PrecompoundModel"); 54 #endif 55 return theReactionProduct; 56 } 57 58 G4double G4PreCompoundNeutron::GetRj(const G4int NumberParticles, const G4int NumberCharged) 46 #include "G4Neutron.hh" 47 48 G4PreCompoundNeutron::G4PreCompoundNeutron() 49 : G4PreCompoundNucleon(G4Neutron::Neutron(), &theNeutronCoulombBarrier) 50 {} 51 52 G4PreCompoundNeutron::~G4PreCompoundNeutron() 53 {} 54 55 G4double G4PreCompoundNeutron::GetRj(G4int nParticles, G4int nCharged) 59 56 { 60 57 G4double rj = 0.0; 61 if(NumberParticles > 0) rj = static_cast<G4double>(NumberParticles - NumberCharged)/ 62 static_cast<G4double>(NumberParticles); 58 if(nParticles > 0) { 59 rj = static_cast<G4double>(nParticles - nCharged)/ 60 static_cast<G4double>(nParticles); 61 } 63 62 return rj; 64 63 } … … 72 71 G4double G4PreCompoundNeutron::CrossSection(const G4double K) 73 72 { 74 ResidualA =GetRestA();75 ResidualZ =GetRestZ();76 theA =GetA();77 theZ =GetZ();78 ResidualAthrd =std::pow(ResidualA,0.33333);79 FragmentA =GetA()+GetRestA();80 FragmentAthrd =std::pow(FragmentA,0.33333);81 82 if (OPTxs==0) return GetOpt0( K);83 else if( OPTxs==1 || OPTxs==2) return GetOpt12( K);84 else if (OPTxs==3 || OPTxs==4) return GetOpt34( K);73 ResidualA = GetRestA(); 74 ResidualZ = GetRestZ(); 75 theA = GetA(); 76 theZ = GetZ(); 77 ResidualAthrd = ResidualA13(); 78 FragmentA = theA + ResidualA; 79 FragmentAthrd = g4pow->Z13(FragmentA); 80 81 if (OPTxs==0) { return GetOpt0( K); } 82 else if( OPTxs==1 || OPTxs==2) { return GetOpt12( K); } 83 else if (OPTxs==3 || OPTxs==4) { return GetOpt34( K); } 85 84 else{ 86 85 std::ostringstream errOs; … … 91 90 } 92 91 93 // *********************** OPT=0 : Dostrovski's cross section *****************************94 95 G4double G4PreCompoundNeutron::GetOpt0(const G4double K)96 {97 98 const G4double r0 = G4PreCompoundParameters::GetAddress()->Getr0();99 // cross section is now given in mb (r0 is in mm) for the sake of consistency100 //with the rest of the options101 return 1.e+25*pi*(r0*ResidualAthrd)*(r0*ResidualAthrd)*GetAlpha()*(1.+GetBeta()/K);102 }103 //104 //-------105 //106 92 G4double G4PreCompoundNeutron::GetAlpha() 107 93 { 108 // return 0.76+2.2/std::pow(GetRestA(),1.0/3.0);109 94 return 0.76+2.2/ResidualAthrd; 110 95 } 111 // 112 //------------ 113 // 96 114 97 G4double G4PreCompoundNeutron::GetBeta() 115 98 { … … 117 100 return (2.12/(ResidualAthrd*ResidualAthrd)-0.05)*MeV/GetAlpha(); 118 101 } 119 // 120 121 //********************* OPT=1,2 : Chatterjee's cross section ************************ 102 103 //********************* OPT=1,2 : Chatterjee's cross section ******************* 122 104 //(fitting to cross section from Bechetti & Greenles OM potential) 123 105 124 G4double G4PreCompoundNeutron::GetOpt12(const G4double K) 125 { 126 106 G4double G4PreCompoundNeutron::GetOpt12(G4double K) 107 { 127 108 G4double Kc=K; 128 109 129 // Pramana (Bechetti & Greenles) for neutrons is chosen130 131 // JMQ xsec is set constat above limit of validity132 if (K>50) Kc=50;110 // Pramana (Bechetti & Greenles) for neutrons is chosen 111 112 // JMQ xsec is set constat above limit of validity 113 if (K > 50*MeV) { Kc = 50*MeV; } 133 114 134 115 G4double landa, landa0, landa1, mu, mu0, mu1,nu, nu0, nu1, nu2,xs; … … 155 136 } 156 137 157 158 138 // *********** OPT=3,4 : Kalbach's cross sections (from PRECO code)************* 159 G4double G4PreCompoundNeutron::GetOpt34(const G4double K) 160 { 161 139 G4double G4PreCompoundNeutron::GetOpt34(G4double K) 140 { 162 141 G4double landa, landa0, landa1, mu, mu0, mu1,nu, nu0, nu1, nu2; 163 142 G4double p, p0, p1, p2; … … 177 156 p2=0.; 178 157 179 180 158 flow = 1.e-18; 181 159 spill= 1.0e+18; 182 160 183 161 // PRECO xs for neutrons is choosen 184 185 162 p0 = -312.; 186 163 p1= 0.; … … 194 171 nu2 = 1280.8; 195 172 196 if (ResidualA < 40 .) signor=0.7+ResidualA*0.0075;197 if (ResidualA > 210 .) signor = 1. + (ResidualA-210.)/250.;173 if (ResidualA < 40) { signor =0.7 + ResidualA*0.0075; } 174 if (ResidualA > 210) { signor = 1. + (ResidualA-210)/250.; } 198 175 landa = landa0/ResidualAthrd + landa1; 199 176 mu = mu0*ResidualAthrd + mu1*ResidualAthrd*ResidualAthrd; … … 201 178 202 179 // JMQ very low energy behaviour corrected (problem for A (apprx.)>60) 203 if (nu < 0.) nu=-nu;180 if (nu < 0.) { nu=-nu; } 204 181 205 182 ec = 0.5; … … 216 193 ecut = 0.; 217 194 cut = a*a - 4.*p*b; 218 if (cut > 0.) ecut = std::sqrt(cut);195 if (cut > 0.) { ecut = std::sqrt(cut); } 219 196 ecut = (ecut-a) / (p+p); 220 197 ecut2 = ecut; 221 if (cut < 0.) ecut2 = ecut - 2.;222 elab = K * FragmentA / ResidualA;198 if (cut < 0.) { ecut2 = ecut - 2.; } 199 elab = K * FragmentA / G4double(ResidualA); 223 200 sig = 0.; 224 201 if (elab <= ec) { //start for E<Ec 225 if (elab > ecut2) sig = (p*elab*elab+a*elab+b) * signor;202 if (elab > ecut2) { sig = (p*elab*elab+a*elab+b) * signor; } 226 203 } //end for E<Ec 227 204 else { //start for E>Ec 228 205 sig = (landa*elab+mu+nu/elab) * signor; 229 206 geom = 0.; 230 if (xnulam < flow || elab < etest) return sig;207 if (xnulam < flow || elab < etest) { return sig; } 231 208 geom = std::sqrt(theA*K); 232 209 geom = 1.23*ResidualAthrd + ra + 4.573/geom; … … 235 212 236 213 } 237 return sig;} 238 239 // ************************** end of cross sections ******************************* 240 241 214 return sig; 215 }
Note: See TracChangeset
for help on using the changeset viewer.