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