Changeset 1347 for trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4TritonEvaporationProbability.cc
- Timestamp:
- Dec 22, 2010, 3:52:27 PM (14 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4TritonEvaporationProbability.cc
r1315 r1347 24 24 // ******************************************************************** 25 25 // 26 //J.M. Quesada (August2008). Based on: 26 // $Id: G4TritonEvaporationProbability.cc,v 1.18 2010/11/17 11:06:03 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-04-ref-00 $ 28 // 29 // J.M. Quesada (August2008). Based on: 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 "G4TritonEvaporationProbability.hh" … … 36 40 G4TritonEvaporationProbability::G4TritonEvaporationProbability() : 37 41 G4EvaporationProbability(3,1,2,&theCoulombBarrier) // A,Z,Gamma,&theCoulombBarrier 38 { 39 40 } 41 42 G4TritonEvaporationProbability::G4TritonEvaporationProbability(const G4TritonEvaporationProbability &) : G4EvaporationProbability() 43 { 44 throw G4HadronicException(__FILE__, __LINE__, "G4TritonEvaporationProbability::copy_constructor meant to not be accessable"); 45 } 46 47 48 49 50 const G4TritonEvaporationProbability & G4TritonEvaporationProbability:: 51 operator=(const G4TritonEvaporationProbability &) 52 { 53 throw G4HadronicException(__FILE__, __LINE__, "G4TritonEvaporationProbability::operator= meant to not be accessable"); 54 return *this; 55 } 56 57 58 G4bool G4TritonEvaporationProbability::operator==(const G4TritonEvaporationProbability &) const 59 { 60 return false; 61 } 62 63 G4bool G4TritonEvaporationProbability::operator!=(const G4TritonEvaporationProbability &) const 64 { 65 return true; 66 } 67 68 G4double G4TritonEvaporationProbability::CalcAlphaParam(const G4Fragment & fragment) 69 { return 1.0 + CCoeficient(static_cast<G4double>(fragment.GetZ()-GetZ()));} 42 {} 43 44 G4TritonEvaporationProbability::~G4TritonEvaporationProbability() 45 {} 46 47 G4double G4TritonEvaporationProbability::CalcAlphaParam(const G4Fragment & fragment) 48 { 49 return 1.0 + CCoeficient(fragment.GetZ_asInt()-GetZ()); 50 } 70 51 71 G4double G4TritonEvaporationProbability::CalcBetaParam(const G4Fragment & ) 72 { return 0.0; } 73 74 G4double G4TritonEvaporationProbability::CCoeficient(const G4double aZ) 75 { 76 // Data comes from 77 // Dostrovsky, Fraenkel and Friedlander 78 // Physical Review, vol 116, num. 3 1959 79 // 80 // const G4int size = 5; 81 // G4double Zlist[5] = { 10.0, 20.0, 30.0, 50.0, 70.0}; 82 // G4double Cp[5] = { 0.50, 0.28, 0.20, 0.15, 0.10}; 83 // C for triton is equal to C for protons divided by 3 84 G4double C = 0.0; 52 G4double G4TritonEvaporationProbability::CalcBetaParam(const G4Fragment & ) 53 { 54 return 0.0; 55 } 56 57 G4double G4TritonEvaporationProbability::CCoeficient(G4int aZ) 58 { 59 // Data comes from 60 // Dostrovsky, Fraenkel and Friedlander 61 // Physical Review, vol 116, num. 3 1959 62 // 63 // const G4int size = 5; 64 // G4double Zlist[5] = { 10.0, 20.0, 30.0, 50.0, 70.0}; 65 // G4double Cp[5] = { 0.50, 0.28, 0.20, 0.15, 0.10}; 66 // C for triton is equal to C for protons divided by 3 67 G4double C = 0.0; 85 68 86 87 88 89 90 69 if (aZ >= 70) { 70 C = 0.10; 71 } else { 72 C = ((((0.15417e-06*aZ) - 0.29875e-04)*aZ + 0.21071e-02)*aZ - 0.66612e-01)*aZ + 0.98375; 73 } 91 74 92 return C/3.0; 93 75 return C/3.0; 94 76 } 95 77 … … 100 82 //OPT=3,4 Kalbach's parameterization 101 83 // 102 G4double G4TritonEvaporationProbability::CrossSection(const G4Fragment & fragment, const G4double K) 84 G4double 85 G4TritonEvaporationProbability::CrossSection(const G4Fragment & fragment, G4double K) 103 86 { 104 87 theA=GetA(); 105 88 theZ=GetZ(); 106 ResidualA=fragment.GetA ()-theA;107 ResidualZ=fragment.GetZ ()-theZ;89 ResidualA=fragment.GetA_asInt()-theA; 90 ResidualZ=fragment.GetZ_asInt()-theZ; 108 91 109 ResidualAthrd=std::pow(ResidualA,0.33333); 110 FragmentA=fragment.GetA(); 111 FragmentAthrd=std::pow(FragmentA,0.33333); 112 113 92 ResidualAthrd=fG4pow->Z13(ResidualA); 93 FragmentA=fragment.GetA_asInt(); 94 FragmentAthrd=fG4pow->Z13(FragmentA); 95 114 96 if (OPTxs==0) {std::ostringstream errOs; 115 97 errOs << "We should'n be here (OPT =0) at evaporation cross section calculation (tritons)!!" … … 128 110 129 111 // 130 //********************* OPT=1,2 : Chatterjee's cross section ***************** *******112 //********************* OPT=1,2 : Chatterjee's cross section ***************** 131 113 //(fitting to cross section from Bechetti & Greenles OM potential) 132 114 133 G4double G4TritonEvaporationProbability::GetOpt12(const G4double K) 134 { 135 115 G4double G4TritonEvaporationProbability::GetOpt12(G4double K) 116 { 136 117 G4double Kc=K; 137 118 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 119 // JMQ xsec is set constat above limit of validity 120 if (K > 50*MeV) { Kc=50*MeV; } 121 122 G4double landa ,mu ,nu ,p , Ec,q,r,ji,xs; 144 123 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 124 G4double p0 = -11.04; 125 G4double p1 = 619.1; 126 G4double p2 = -2147.; 127 G4double landa0 = -0.0426; 128 G4double landa1 = -10.33; 129 G4double mu0 = 601.9; 130 G4double mu1 = 0.37; 131 G4double nu0 = 583.0; 132 G4double nu1 = -546.2; 133 G4double nu2 = 1.718; 134 G4double delta=1.2; 135 136 Ec = 1.44*theZ*ResidualZ/(1.5*ResidualAthrd+delta); 137 p = p0 + p1/Ec + p2/(Ec*Ec); 138 landa = landa0*ResidualA + landa1; 139 140 G4double resmu1 = fG4pow->powZ(ResidualA,mu1); 141 mu = mu0*resmu1; 142 nu = resmu1*(nu0 + nu1*Ec + nu2*(Ec*Ec)); 143 q = landa - nu/(Ec*Ec) - 2*p*Ec; 144 r = mu + 2*nu/Ec + p*(Ec*Ec); 145 146 ji=std::max(Kc,Ec); 147 if(Kc < Ec) { xs = p*Kc*Kc + q*Kc + r;} 148 else {xs = p*(Kc - ji)*(Kc - ji) + landa*Kc + mu + nu*(2 - Kc/ji)/ji ;} 149 150 if (xs <0.0) {xs=0.0;} 151 152 return xs; 173 153 } 174 154 175 155 // *********** OPT=3,4 : Kalbach's cross sections (from PRECO code)************* 176 G4double G4TritonEvaporationProbability::GetOpt34( constG4double K)156 G4double G4TritonEvaporationProbability::GetOpt34(G4double K) 177 157 // ** t from o.m. of hafele, flynn et al 178 158 { 179 180 159 G4double landa, mu, nu, p , signor(1.),sig; 181 160 G4double ec,ecsq,xnulam,etest(0.),a; 182 161 G4double b,ecut,cut,ecut2,geom,elab; 183 184 162 185 163 G4double flow = 1.e-18; 186 164 G4double spill= 1.e+18; … … 205 183 p = p0 + p1/ec + p2/ecsq; 206 184 landa = landa0*ResidualA + landa1; 207 a = std::pow(ResidualA,mu1);185 a = fG4pow->powZ(ResidualA,mu1); 208 186 mu = mu0 * a; 209 187 nu = a* (nu0+nu1*ec+nu2*ecsq); 210 188 xnulam = nu / landa; 211 if (xnulam > spill) xnulam=0.;212 if (xnulam >= flow) etest = 1.2 *std::sqrt(xnulam);213 189 if (xnulam > spill) { xnulam=0.; } 190 if (xnulam >= flow) { etest = 1.2 *std::sqrt(xnulam); } 191 214 192 a = -2.*p*ec + landa - nu/ecsq; 215 193 b = p*ecsq + mu + 2.*nu/ec; 216 194 ecut = 0.; 217 195 cut = a*a - 4.*p*b; 218 if (cut > 0.) ecut = std::sqrt(cut);196 if (cut > 0.) { ecut = std::sqrt(cut); } 219 197 ecut = (ecut-a) / (p+p); 220 198 ecut2 = ecut; 221 //JMQ 290310 for avoiding unphysical increase below minimum (at ecut) 222 //ecut<0 means that there is no cut with energy axis, i.e. xs is set to 0 bellow minimum 223 // if (cut < 0.) ecut2 = ecut - 2.; 224 if (cut < 0.) ecut2 = ecut; 225 elab = K * FragmentA / ResidualA; 199 //JMQ 290310 for avoiding unphysical increase below minimum (at ecut) 200 // ecut<0 means that there is no cut with energy axis, i.e. xs is set 201 // to 0 bellow minimum 202 // if (cut < 0.) ecut2 = ecut - 2.; 203 if (cut < 0.) { ecut2 = ecut; } 204 elab = K * FragmentA / G4double(ResidualA); 226 205 sig = 0.; 227 206 228 207 if (elab <= ec) { //start for E<Ec 229 if (elab > ecut2) sig = (p*elab*elab+a*elab+b) * signor;208 if (elab > ecut2) { sig = (p*elab*elab+a*elab+b) * signor; } 230 209 } //end for E<Ec 231 else { //start for E>Ec 210 else { //start for E>Ec 232 211 sig = (landa*elab+mu+nu/elab) * signor; 233 212 geom = 0.; 234 if (xnulam < flow || elab < etest) return sig;213 if (xnulam < flow || elab < etest) { return sig; } 235 214 geom = std::sqrt(theA*K); 236 215 geom = 1.23*ResidualAthrd + ra + 4.573/geom; … … 239 218 } //end for E>Ec 240 219 return sig; 241 242 } 243 244 // ************************** end of cross sections ******************************* 245 220 } 221
Note: See TracChangeset
for help on using the changeset viewer.