- Timestamp:
- Dec 22, 2010, 3:52:27 PM (14 years ago)
- Location:
- trunk/source/processes/hadronic/models/de_excitation/evaporation/src
- Files:
-
- 19 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4AlphaEvaporationChannel.cc
r1340 r1347 25 25 // 26 26 // 27 // $Id: G4AlphaEvaporationChannel.cc,v 1. 4 2006/06/29 20:10:17 gunterExp $28 // GEANT4 tag $Name: geant4-09-0 3-ref-09$27 // $Id: G4AlphaEvaporationChannel.cc,v 1.5 2010/11/17 12:14:59 vnivanch Exp $ 28 // GEANT4 tag $Name: geant4-09-04-ref-00 $ 29 29 // 30 30 // Hadronic Process: Nuclear De-excitations 31 31 // by V. Lara (Nov. 1999) 32 32 // 33 // 17-11-2010 V.Ivanchenko moved constructor and destructor to source and cleanup 33 34 34 35 #include "G4AlphaEvaporationChannel.hh" 35 36 37 G4AlphaEvaporationChannel::G4AlphaEvaporationChannel() 38 : G4EvaporationChannel(4,2,"alpha",&theEvaporationProbability,&theCoulombBarrier) 39 {} 36 40 37 const G4AlphaEvaporationChannel & G4AlphaEvaporationChannel::operator=(const G4AlphaEvaporationChannel & ) 38 { 39 throw G4HadronicException(__FILE__, __LINE__, "G4AlphaEvaporationChannel::operator= meant to not be accessable"); 40 return *this; 41 } 41 G4AlphaEvaporationChannel::~G4AlphaEvaporationChannel() 42 {} 42 43 43 G4AlphaEvaporationChannel::G4AlphaEvaporationChannel(const G4AlphaEvaporationChannel & ) : G4EvaporationChannel()44 {45 throw G4HadronicException(__FILE__, __LINE__, "G4AlphaEvaporationChannel::CopyConstructor meant to not be accessable");46 }47 48 G4bool G4AlphaEvaporationChannel::operator==(const G4AlphaEvaporationChannel &right ) const49 {50 return (this == (G4AlphaEvaporationChannel *) &right);51 // return false;52 }53 54 G4bool G4AlphaEvaporationChannel::operator!=(const G4AlphaEvaporationChannel &right ) const55 {56 return (this != (G4AlphaEvaporationChannel *) &right);57 // return true;58 }59 -
trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4AlphaEvaporationProbability.cc
r1315 r1347 24 24 // ******************************************************************** 25 25 // 26 //J.M. Quesada (August2008). Based on: 26 // $Id: G4AlphaEvaporationProbability.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 "G4AlphaEvaporationProbability.hh" … … 36 40 G4AlphaEvaporationProbability::G4AlphaEvaporationProbability() : 37 41 G4EvaporationProbability(4,2,1,&theCoulombBarrier) // A,Z,Gamma,&theCoumlombBarrier 38 { 42 {} 43 44 G4AlphaEvaporationProbability::~G4AlphaEvaporationProbability() 45 {} 46 47 G4double G4AlphaEvaporationProbability::CalcAlphaParam(const G4Fragment & fragment) 48 { return 1.0 + CCoeficient(fragment.GetZ_asInt()-GetZ());} 39 49 40 } 41 42 43 G4AlphaEvaporationProbability::G4AlphaEvaporationProbability(const G4AlphaEvaporationProbability &): G4EvaporationProbability() 44 { 45 throw G4HadronicException(__FILE__, __LINE__, "G4AlphaEvaporationProbability::copy_constructor meant to not be accessable"); 46 } 47 48 49 50 51 const G4AlphaEvaporationProbability & G4AlphaEvaporationProbability:: 52 operator=(const G4AlphaEvaporationProbability &) 53 { 54 throw G4HadronicException(__FILE__, __LINE__, "G4AlphaEvaporationProbability::operator= meant to not be accessable"); 55 return *this; 56 } 57 58 59 G4bool G4AlphaEvaporationProbability::operator==(const G4AlphaEvaporationProbability &) const 60 { 61 return false; 62 } 63 64 G4bool G4AlphaEvaporationProbability::operator!=(const G4AlphaEvaporationProbability &) const 65 { 66 return true; 67 } 68 69 70 G4double G4AlphaEvaporationProbability::CalcAlphaParam(const G4Fragment & fragment) 71 { return 1.0 + CCoeficient(static_cast<G4double>(fragment.GetZ()-GetZ()));} 50 G4double G4AlphaEvaporationProbability::CalcBetaParam(const G4Fragment &) 51 { return 0.0; } 52 53 G4double G4AlphaEvaporationProbability::CCoeficient(G4int aZ) 54 { 55 // Data comes from 56 // Dostrovsky, Fraenkel and Friedlander 57 // Physical Review, vol 116, num. 3 1959 58 // 59 // const G4int size = 5; 60 // G4double Zlist[5] = { 10.0, 20.0, 30.0, 50.0, 70.0}; 61 // G4double Calpha[5] = { 0.10, 0.10, 0.10, 0.08, 0.06}; 62 G4double C = 0.0; 72 63 73 G4double G4AlphaEvaporationProbability::CalcBetaParam(const G4Fragment &) 74 { return 0.0; } 75 76 G4double G4AlphaEvaporationProbability::CCoeficient(const G4double aZ) 77 { 78 // Data comes from 79 // Dostrovsky, Fraenkel and Friedlander 80 // Physical Review, vol 116, num. 3 1959 81 // 82 // const G4int size = 5; 83 // G4double Zlist[5] = { 10.0, 20.0, 30.0, 50.0, 70.0}; 84 // G4double Calpha[5] = { 0.10, 0.10, 0.10, 0.08, 0.06}; 85 G4double C = 0.0; 86 87 88 if (aZ <= 30) { 89 C = 0.10; 90 } else if (aZ <= 50) { 91 C = 0.1 + -((aZ-50.)/20.)*0.02; 92 } else if (aZ < 70) { 93 C = 0.08 + -((aZ-70.)/20.)*0.02; 94 } else { 95 C = 0.06; 96 } 97 return C; 64 if (aZ <= 30) 65 { 66 C = 0.10; 67 } 68 else if (aZ <= 50) 69 { 70 C = 0.1 - (aZ-30)*0.001; 71 } 72 else if (aZ < 70) 73 { 74 C = 0.08 - (aZ-50)*0.001; 75 } 76 else 77 { 78 C = 0.06; 79 } 80 return C; 98 81 } 99 82 … … 104 87 //OPT=3,4 Kalbach's parameterization 105 88 // 106 G4double G4AlphaEvaporationProbability::CrossSection(const G4Fragment & fragment,107 constG4double K)89 G4double 90 G4AlphaEvaporationProbability::CrossSection(const G4Fragment & fragment, G4double K) 108 91 { 109 92 theA=GetA(); 110 93 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 94 ResidualA=fragment.GetA_asInt()-theA; 95 ResidualZ=fragment.GetZ_asInt()-theZ; 96 97 ResidualAthrd=fG4pow->Z13(ResidualA); 98 FragmentA=fragment.GetA_asInt(); 99 FragmentAthrd=fG4pow->Z13(FragmentA); 118 100 119 101 if (OPTxs==0) {std::ostringstream errOs; … … 133 115 } 134 116 135 136 // 137 //********************* OPT=1,2 : Chatterjee's cross section ************************ 117 // 118 //********************* OPT=1,2 : Chatterjee's cross section ******************** 138 119 //(fitting to cross section from Bechetti & Greenles OM potential) 139 120 140 G4double G4AlphaEvaporationProbability::GetOpt12(const G4double K) 141 { 142 // c ** alpha from huizenga and igo 121 G4double G4AlphaEvaporationProbability::GetOpt12(G4double K) 122 { 143 123 G4double Kc=K; 144 145 // JMQ xsec is set consta t above limit of validity146 if (K >50) Kc=50;147 124 125 // JMQ xsec is set constant above limit of validity 126 if (K > 50*MeV) { Kc = 50*MeV; } 127 148 128 G4double landa ,mu ,nu ,p , Ec,q,r,ji,xs; 149 129 150 130 G4double p0 = 10.95; 151 131 G4double p1 = -85.2; … … 160 140 G4double delta=1.2; 161 141 162 163 142 Ec = 1.44*theZ*ResidualZ/(1.5*ResidualAthrd+delta); 164 143 p = p0 + p1/Ec + p2/(Ec*Ec); 165 144 landa = landa0*ResidualA + landa1; 166 mu = mu0*std::pow(ResidualA,mu1); 167 nu = std::pow(ResidualA,mu1)*(nu0 + nu1*Ec + nu2*(Ec*Ec)); 145 G4double resmu1 = fG4pow->powZ(ResidualA,mu1); 146 mu = mu0*resmu1; 147 nu = resmu1*(nu0 + nu1*Ec + nu2*(Ec*Ec)); 168 148 q = landa - nu/(Ec*Ec) - 2*p*Ec; 169 149 r = mu + 2*nu/Ec + p*(Ec*Ec); 170 150 171 151 ji=std::max(Kc,Ec); 172 152 if(Kc < Ec) { xs = p*Kc*Kc + q*Kc + r;} 173 153 else {xs = p*(Kc - ji)*(Kc - ji) + landa*Kc + mu + nu*(2 - Kc/ji)/ji ;} 174 154 175 155 if (xs <0.0) {xs=0.0;} 176 156 177 157 return xs; 178 179 158 } 180 159 181 160 // *********** OPT=3,4 : Kalbach's cross sections (from PRECO code)************* 182 G4double G4AlphaEvaporationProbability::GetOpt34( constG4double K)161 G4double G4AlphaEvaporationProbability::GetOpt34(G4double K) 183 162 // c ** alpha from huizenga and igo 184 163 { 185 186 164 G4double landa, mu, nu, p , signor(1.),sig; 187 165 G4double ec,ecsq,xnulam,etest(0.),a; 188 166 G4double b,ecut,cut,ecut2,geom,elab; 189 167 190 191 168 G4double flow = 1.e-18; 192 169 G4double spill= 1.e+18; 193 170 194 G4double 171 G4double p0 = 10.95; 195 172 G4double p1 = -85.2; 196 173 G4double p2 = 1146.; … … 204 181 205 182 G4double ra=1.20; 206 183 207 184 //JMQ 13/02/09 increase of reduced radius to lower the barrier 208 185 // ec = 1.44 * theZ * ResidualZ / (1.5*ResidualAthrd+ra); … … 211 188 p = p0 + p1/ec + p2/ecsq; 212 189 landa = landa0*ResidualA + landa1; 213 a = std::pow(ResidualA,mu1);190 a = fG4pow->powZ(ResidualA,mu1); 214 191 mu = mu0 * a; 215 192 nu = a* (nu0+nu1*ec+nu2*ecsq); 216 193 xnulam = nu / landa; 217 if (xnulam > spill) xnulam=0.;218 if (xnulam >= flow) etest = 1.2 *std::sqrt(xnulam);219 194 if (xnulam > spill) { xnulam=0.; } 195 if (xnulam >= flow) { etest = 1.2 *std::sqrt(xnulam); } 196 220 197 a = -2.*p*ec + landa - nu/ecsq; 221 198 b = p*ecsq + mu + 2.*nu/ec; 222 199 ecut = 0.; 223 200 cut = a*a - 4.*p*b; 224 if (cut > 0.) ecut = std::sqrt(cut);201 if (cut > 0.) { ecut = std::sqrt(cut); } 225 202 ecut = (ecut-a) / (p+p); 226 203 ecut2 = ecut; 227 //JMQ 290310 for avoiding unphysical increase below minimum (at ecut) 228 //ecut<0 means that there is no cut with energy axis, i.e. xs is set to 0 bellow minimum 229 // if (cut < 0.) ecut2 = ecut - 2.; 230 if (cut < 0.) ecut2 = ecut; 231 elab = K * FragmentA / ResidualA; 204 //JMQ 290310 for avoiding unphysical increase below minimum (at ecut) 205 // ecut<0 means that there is no cut with energy axis, i.e. xs is set 206 // to 0 bellow minimum 207 // if (cut < 0.) ecut2 = ecut - 2.; 208 if (cut < 0.) { ecut2 = ecut; } 209 elab = K * FragmentA / G4double(ResidualA); 232 210 sig = 0.; 233 211 234 212 if (elab <= ec) { //start for E<Ec 235 if (elab > ecut2) sig = (p*elab*elab+a*elab+b) * signor;213 if (elab > ecut2) { sig = (p*elab*elab+a*elab+b) * signor; } 236 214 } //end for E<Ec 237 215 else { //start for E>Ec 238 216 sig = (landa*elab+mu+nu/elab) * signor; 239 217 geom = 0.; 240 if (xnulam < flow || elab < etest) return sig;218 if (xnulam < flow || elab < etest) { return sig; } 241 219 geom = std::sqrt(theA*K); 242 220 geom = 1.23*ResidualAthrd + ra + 4.573/geom; … … 245 223 } //end for E>Ec 246 224 return sig; 247 248 } 249 250 // ************************** end of cross sections ******************************* 225 } 226 -
trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4DeuteronEvaporationChannel.cc
r1340 r1347 25 25 // 26 26 // 27 // $Id: G4DeuteronEvaporationChannel.cc,v 1. 4 2006/06/29 20:10:21 gunterExp $28 // GEANT4 tag $Name: geant4-09-0 3-ref-09$27 // $Id: G4DeuteronEvaporationChannel.cc,v 1.5 2010/11/17 12:14:59 vnivanch Exp $ 28 // GEANT4 tag $Name: geant4-09-04-ref-00 $ 29 29 // 30 30 // Hadronic Process: Nuclear De-excitations 31 31 // by V. Lara (Nov. 1999) 32 32 // 33 // 17-11-2010 V.Ivanchenko moved constructor and destructor to source and cleanup 33 34 34 35 #include "G4DeuteronEvaporationChannel.hh" 35 36 37 G4DeuteronEvaporationChannel::G4DeuteronEvaporationChannel() 38 : G4EvaporationChannel(2,1,"deuteron",&theEvaporationProbability,&theCoulombBarrier) 39 {} 36 40 37 const G4DeuteronEvaporationChannel & G4DeuteronEvaporationChannel::operator=(const G4DeuteronEvaporationChannel & ) 38 { 39 throw G4HadronicException(__FILE__, __LINE__, "G4DeuteronEvaporationChannel::operator= meant to not be accessable"); 40 return *this; 41 } 41 G4DeuteronEvaporationChannel::~G4DeuteronEvaporationChannel() 42 {} 42 43 43 G4DeuteronEvaporationChannel::G4DeuteronEvaporationChannel(const G4DeuteronEvaporationChannel & ) : G4EvaporationChannel()44 {45 throw G4HadronicException(__FILE__, __LINE__, "G4DeuteronEvaporationChannel::CopyConstructor meant to not be accessable");46 }47 48 G4bool G4DeuteronEvaporationChannel::operator==(const G4DeuteronEvaporationChannel & right) const49 {50 return (this == (G4DeuteronEvaporationChannel *) &right);51 // return false;52 }53 54 G4bool G4DeuteronEvaporationChannel::operator!=(const G4DeuteronEvaporationChannel & right) const55 {56 return (this != (G4DeuteronEvaporationChannel *) &right);57 // return true;58 } -
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 -
trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4Evaporation.cc
r1340 r1347 25 25 // 26 26 // 27 // $Id: G4Evaporation.cc,v 1.2 5 2010/06/09 11:56:47vnivanch Exp $28 // GEANT4 tag $Name: geant4-09-0 3-ref-09$27 // $Id: G4Evaporation.cc,v 1.26 2010/11/23 18:10:10 vnivanch Exp $ 28 // GEANT4 tag $Name: geant4-09-04-ref-00 $ 29 29 // 30 30 // Hadronic Process: Nuclear De-excitations … … 142 142 G4int Z = theResidualNucleus->GetZ_asInt(); 143 143 G4double abun = nist->GetIsotopeAbundance(Z, A); 144 /*145 G4cout << "### G4Evaporation::BreakItUp step " << ia << " Z= " << Z146 << " A= " << A << " Eex(MeV)= "147 << theResidualNucleus->GetExcitationEnergy()148 << " aban= " << abun << G4endl;149 */144 145 // G4cout << "### G4Evaporation::BreakItUp step " << ia << " Z= " << Z 146 // << " A= " << A << " Eex(MeV)= " 147 // << theResidualNucleus->GetExcitationEnergy() 148 // << " aban= " << abun << G4endl; 149 150 150 if(theResidualNucleus->GetExcitationEnergy() <= minExcitation && 151 151 (abun > 0.0)) … … 184 184 // do evaporation chain and reset total probability 185 185 if(0.0 < totprob && probabilities[0] == totprob) { 186 //G4cout << "Start gamma evaporation" << G4endl; 186 187 theTempResult = (*theChannels)[0]->BreakUpFragment(theResidualNucleus); 187 188 if(theTempResult) { … … 228 229 // single photon evaporation, primary pointer is kept 229 230 if(0 == i) { 231 //G4cout << "Single gamma" << G4endl; 230 232 G4Fragment* gamma = (*theChannels)[0]->EmittedFragment(theResidualNucleus); 231 233 if(gamma) { theResult->push_back(gamma); } … … 233 235 // fission, return results to the main loop if fission is succesful 234 236 } else if(1 == i) { 237 //G4cout << "Fission" << G4endl; 235 238 theTempResult = (*theChannels)[1]->BreakUp(*theResidualNucleus); 236 239 if(theTempResult) { … … 248 251 // other channels 249 252 } else { 253 //G4cout << "Channel # " << i << G4endl; 250 254 theTempResult = (*theChannels)[i]->BreakUp(*theResidualNucleus); 251 255 if(theTempResult) { -
trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4EvaporationChannel.cc
r962 r1347 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4EvaporationChannel.cc,v 1.19 2010/11/24 15:30:49 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-04-ref-00 $ 26 28 // 27 29 //J.M. Quesada (August2008). Based on: … … 30 32 // by V. Lara (Oct 1998) 31 33 // 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 34 // Modified: 35 // 03-09-2008 J.M. Quesada for external choice of inverse cross section option 36 // 06-09-2008 J.M. Quesada Also external choices have been added for superimposed 37 // Coulomb barrier (if useSICB is set true, by default is false) 38 // 17-11-2010 V.Ivanchenko in constructor replace G4VEmissionProbability by 39 // G4EvaporationProbability and do not new and delete probability 40 // object at each call; use G4Pow 37 41 38 42 #include "G4EvaporationChannel.hh" 39 43 #include "G4PairingCorrection.hh" 40 41 42 43 G4EvaporationChannel::G4EvaporationChannel(const G4int anA, const G4int aZ, const G4String & aName, 44 G4VEmissionProbability * aEmissionStrategy, 45 G4VCoulombBarrier * aCoulombBarrier): 44 #include "G4NucleiProperties.hh" 45 #include "G4Pow.hh" 46 #include "G4EvaporationLevelDensityParameter.hh" 47 #include "Randomize.hh" 48 #include "G4Alpha.hh" 49 50 G4EvaporationChannel::G4EvaporationChannel(G4int anA, G4int aZ, 51 const G4String & aName, 52 G4EvaporationProbability* aEmissionStrategy, 53 G4VCoulombBarrier* aCoulombBarrier): 46 54 G4VEvaporationChannel(aName), 47 55 theA(anA), … … 52 60 MaximalKineticEnergy(-1000.0) 53 61 { 54 theLevelDensityPtr = new G4EvaporationLevelDensityParameter; 55 // MyOwnLevelDensity = true; 62 ResidualA = 0; 63 ResidualZ = 0; 64 ResidualMass = CoulombBarrier=0.0; 65 EvaporatedMass = G4NucleiProperties::GetNuclearMass(theA, theZ); 66 theLevelDensityPtr = new G4EvaporationLevelDensityParameter; 67 } 68 69 G4EvaporationChannel::G4EvaporationChannel(): 70 G4VEvaporationChannel(""), 71 theA(0), 72 theZ(0), 73 theEvaporationProbabilityPtr(0), 74 theCoulombBarrierPtr(0), 75 EmissionProbability(0.0), 76 MaximalKineticEnergy(-1000.0) 77 { 78 ResidualA = 0; 79 ResidualZ = 0; 80 EvaporatedMass = ResidualMass = CoulombBarrier = 0.0; 81 theLevelDensityPtr = new G4EvaporationLevelDensityParameter; 56 82 } 57 83 … … 60 86 delete theLevelDensityPtr; 61 87 } 62 63 G4EvaporationChannel::G4EvaporationChannel(const G4EvaporationChannel & ) : G4VEvaporationChannel()64 {65 throw G4HadronicException(__FILE__, __LINE__, "G4EvaporationChannel::copy_costructor meant to not be accessable");66 }67 68 const G4EvaporationChannel & G4EvaporationChannel::operator=(const G4EvaporationChannel & )69 {70 throw G4HadronicException(__FILE__, __LINE__, "G4EvaporationChannel::operator= meant to not be accessable");71 return *this;72 }73 74 G4bool G4EvaporationChannel::operator==(const G4EvaporationChannel & right) const75 {76 return (this == (G4EvaporationChannel *) &right);77 // return false;78 }79 80 G4bool G4EvaporationChannel::operator!=(const G4EvaporationChannel & right) const81 {82 return (this != (G4EvaporationChannel *) &right);83 // return true;84 }85 86 //void G4EvaporationChannel::SetLevelDensityParameter(G4VLevelDensityParameter * aLevelDensity)87 // {88 // if (MyOwnLevelDensity) delete theLevelDensityPtr;89 // theLevelDensityPtr = aLevelDensity;90 // MyOwnLevelDensity = false;91 // }92 88 93 89 void G4EvaporationChannel::Initialize(const G4Fragment & fragment) … … 98 94 theEvaporationProbabilityPtr->UseSICB(useSICB); 99 95 100 101 G4int FragmentA = static_cast<G4int>(fragment.GetA()); 102 G4int FragmentZ = static_cast<G4int>(fragment.GetZ()); 96 G4int FragmentA = fragment.GetA_asInt(); 97 G4int FragmentZ = fragment.GetZ_asInt(); 103 98 ResidualA = FragmentA - theA; 104 99 ResidualZ = FragmentZ - theZ; 100 //G4cout << "G4EvaporationChannel::Initialize Z= " << theZ << " A= " << theA 101 // << " FragZ= " << FragmentZ << " FragA= " << FragmentA << G4endl; 105 102 106 103 //Effective excitation energy … … 108 105 G4PairingCorrection::GetInstance()->GetPairingCorrection(FragmentA,FragmentZ); 109 106 110 111 107 // Only channels which are physically allowed are taken into account 112 108 if (ResidualA <= 0 || ResidualZ <= 0 || ResidualA < ResidualZ || 113 109 (ResidualA == ResidualZ && ResidualA > 1) || ExEnergy <= 0.0) { 114 CoulombBarrier =0.0;110 CoulombBarrier = ResidualMass = 0.0; 115 111 MaximalKineticEnergy = -1000.0*MeV; 116 112 EmissionProbability = 0.0; 117 113 } else { 114 ResidualMass = G4NucleiProperties::GetNuclearMass(ResidualA, ResidualZ); 115 G4double FragmentMass = fragment.GetGroundStateMass(); 118 116 CoulombBarrier = theCoulombBarrierPtr->GetCoulombBarrier(ResidualA,ResidualZ,ExEnergy); 119 117 // Maximal Kinetic Energy 120 MaximalKineticEnergy = CalcMaximalKineticEnergy 121 (G4ParticleTable::GetParticleTable()->122 GetIonTable()->GetNucleusMass(FragmentZ,FragmentA)+ExEnergy);118 MaximalKineticEnergy = CalcMaximalKineticEnergy(FragmentMass + ExEnergy); 119 //MaximalKineticEnergy = ExEnergy + fragment.GetGroundStateMass() 120 // - EvaporatedMass - ResidualMass; 123 121 124 122 // Emission probability … … 131 129 limit= CoulombBarrier; 132 130 else limit=0.; 131 limit = std::max(limit, FragmentMass - ResidualMass - EvaporatedMass); 133 132 134 133 // The threshold for charged particle emission must be set to 0 if Coulomb 135 134 //cutoff is included in the cross sections 136 //if (MaximalKineticEnergy <= 0.0) EmissionProbability = 0.0;137 135 if (MaximalKineticEnergy <= limit) EmissionProbability = 0.0; 138 136 else { … … 142 140 } 143 141 } 142 //G4cout << "G4EvaporationChannel:: probability= " << EmissionProbability << G4endl; 144 143 145 144 return; … … 148 147 G4FragmentVector * G4EvaporationChannel::BreakUp(const G4Fragment & theNucleus) 149 148 { 150 G4double EvaporatedKineticEnergy=GetKineticEnergy(theNucleus); 151 152 G4double EvaporatedMass = G4ParticleTable::GetParticleTable()->GetIonTable()-> 153 GetIonMass(theZ,theA); 154 G4double EvaporatedEnergy = EvaporatedKineticEnergy + EvaporatedMass; 155 149 /* 150 G4double Ecm = GetKineticEnergy(theNucleus) + ResidualMass + EvaporatedMass; 151 152 G4double EvaporatedEnergy = 153 ((Ecm-ResidualMass)*(Ecm+ResidualMass) + EvaporatedMass*EvaporatedMass)/(2*Ecm); 154 */ 155 G4double EvaporatedEnergy = GetKineticEnergy(theNucleus) + EvaporatedMass; 156 156 157 G4ThreeVector momentum(IsotropicVector 157 (std::sqrt(EvaporatedKineticEnergy* 158 (EvaporatedKineticEnergy+2.0*EvaporatedMass)))); 159 160 momentum.rotateUz(theNucleus.GetMomentum().vect().unit()); 158 (std::sqrt((EvaporatedEnergy - EvaporatedMass)* 159 (EvaporatedEnergy + EvaporatedMass)))); 161 160 162 161 G4LorentzVector EvaporatedMomentum(momentum,EvaporatedEnergy); 163 EvaporatedMomentum.boost(theNucleus.GetMomentum().boostVector()); 162 G4LorentzVector ResidualMomentum = theNucleus.GetMomentum(); 163 EvaporatedMomentum.boost(ResidualMomentum.boostVector()); 164 164 165 165 G4Fragment * EvaporatedFragment = new G4Fragment(theA,theZ,EvaporatedMomentum); 166 #ifdef PRECOMPOUND_TEST 167 EvaporatedFragment->SetCreatorModel(G4String("G4Evaporation")); 168 #endif 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 182 #ifdef PRECOMPOUND_TEST 183 ResidualFragment->SetCreatorModel(G4String("ResidualNucleus")); 184 #endif 166 ResidualMomentum -= EvaporatedMomentum; 167 168 G4Fragment * ResidualFragment = new G4Fragment(ResidualA, ResidualZ, ResidualMomentum); 169 185 170 G4FragmentVector * theResult = new G4FragmentVector; 186 171 … … 211 196 ///////////////////////////////////////// 212 197 // Calculates the maximal kinetic energy that can be carried by fragment. 213 G4double G4EvaporationChannel::CalcMaximalKineticEnergy(const G4double NucleusTotalE) 214 { 215 G4double ResidualMass = G4ParticleTable::GetParticleTable()-> 216 GetIonTable()->GetNucleusMass( ResidualZ, ResidualA ); 217 G4double EvaporatedMass = G4ParticleTable::GetParticleTable()-> 218 GetIonTable()->GetNucleusMass( theZ, theA ); 219 198 G4double G4EvaporationChannel::CalcMaximalKineticEnergy(G4double NucleusTotalE) 199 { 220 200 // This is the "true" assimptotic kinetic energy (from energy conservation) 221 G4double Tmax = (NucleusTotalE*NucleusTotalE + EvaporatedMass*EvaporatedMass - 222 ResidualMass*ResidualMass)/(2.0*NucleusTotalE) - EvaporatedMass; 201 G4double Tmax = 202 ((NucleusTotalE-ResidualMass)*(NucleusTotalE+ResidualMass) + EvaporatedMass*EvaporatedMass) 203 /(2.0*NucleusTotalE) - EvaporatedMass; 223 204 224 205 //JMQ (13-09-08) bug fixed: in the original version the Tmax is calculated … … 245 226 246 227 247 if (MaximalKineticEnergy < 0.0) 228 if (MaximalKineticEnergy < 0.0) { 248 229 throw G4HadronicException(__FILE__, __LINE__, 249 230 "G4EvaporationChannel::CalcKineticEnergy: maximal kinetic at the Coulomb barrier is less than 0"); 250 231 } 251 232 G4double Rb = 4.0*theLevelDensityPtr-> 252 233 LevelDensityParameter(ResidualA+theA,ResidualZ+theZ,MaximalKineticEnergy)* … … 263 244 G4double Q2 = 1.0; 264 245 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.));246 G4double Beta = (2.12/G4Pow::GetInstance()->Z23(ResidualA) - 0.05)*MeV/ 247 (0.76 + 2.2/G4Pow::GetInstance()->Z13(ResidualA)); 267 248 Q1 = 1.0 + Beta/(MaximalKineticEnergy); 268 249 Q2 = Q1*std::sqrt(Q1); … … 277 258 } else if (OPTxs==1 || OPTxs==2 || OPTxs==3 || OPTxs==4) { 278 259 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.;260 // Coulomb barrier is just included in the cross sections 261 G4double V = 0; 262 if(useSICB) { V= CoulombBarrier; } 263 264 V = std::max(V, aFragment.GetGroundStateMass()-EvaporatedMass-ResidualMass); 284 265 285 266 G4double Tmax=MaximalKineticEnergy; 286 267 G4double T(0.0); 287 268 G4double NormalizedProbability(1.0); 288 269 270 // VI: This is very ineffective - create new objects at each call to the method 271 /* 289 272 // A pointer is created in order to access the distribution function. 290 G4EvaporationProbability * G4EPtemp ;273 G4EvaporationProbability * G4EPtemp = 0; 291 274 292 275 if (theA==1 && theZ==0) G4EPtemp=new G4NeutronEvaporationProbability(); … … 305 288 G4EPtemp->SetOPTxs(OPTxs); 306 289 G4EPtemp->UseSICB(useSICB); 307 290 */ 291 292 // use local pointer and not create a new one 308 293 do 309 294 { 310 295 T=V+G4UniformRand()*(Tmax-V); 311 NormalizedProbability=G4EPtemp->ProbabilityDistributionFunction(aFragment,T)/ 312 (this->GetEmissionProbability()); 296 NormalizedProbability = 297 theEvaporationProbabilityPtr->ProbabilityDistributionFunction(aFragment,T)/ 298 GetEmissionProbability(); 313 299 314 300 } 315 301 while (G4UniformRand() > NormalizedProbability); 316 delete G4EPtemp;317 return T;302 // delete G4EPtemp; 303 return T; 318 304 } else{ 319 305 std::ostringstream errOs; … … 323 309 } 324 310 325 G4ThreeVector G4EvaporationChannel::IsotropicVector( constG4double Magnitude)311 G4ThreeVector G4EvaporationChannel::IsotropicVector(G4double Magnitude) 326 312 // Samples a isotropic random vectorwith a magnitud given by Magnitude. 327 313 // By default Magnitude = 1.0 328 314 { 329 330 331 332 333 334 335 336 315 G4double CosTheta = 1.0 - 2.0*G4UniformRand(); 316 G4double SinTheta = std::sqrt(1.0 - CosTheta*CosTheta); 317 G4double Phi = twopi*G4UniformRand(); 318 G4ThreeVector Vector(Magnitude*std::cos(Phi)*SinTheta, 319 Magnitude*std::sin(Phi)*SinTheta, 320 Magnitude*CosTheta); 321 return Vector; 322 } 337 323 338 324 -
trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4EvaporationDefaultGEMFactory.cc
r1340 r1347 26 26 // 27 27 // $Id: G4EvaporationDefaultGEMFactory.cc,v 1.2 2010/04/27 11:43:16 vnivanch Exp $ 28 // GEANT4 tag $Name: geant4-09-0 3-ref-09$28 // GEANT4 tag $Name: geant4-09-04-ref-00 $ 29 29 // 30 30 // Hadronic Process: Nuclear De-excitations -
trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4EvaporationFactory.cc
r1340 r1347 26 26 // 27 27 // $Id: G4EvaporationFactory.cc,v 1.5 2010/04/27 11:43:16 vnivanch Exp $ 28 // GEANT4 tag $Name: geant4-09-0 3-ref-09$28 // GEANT4 tag $Name: geant4-09-04-ref-00 $ 29 29 // 30 30 // Hadronic Process: Nuclear De-excitations -
trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4EvaporationProbability.cc
r1315 r1347 44 44 #include "G4IonTable.hh" 45 45 46 47 G4EvaporationProbability::G4EvaporationProbability(const G4EvaporationProbability &) : G4VEmissionProbability() 46 G4EvaporationProbability::G4EvaporationProbability(G4int anA, G4int aZ, 47 G4double aGamma, 48 G4VCoulombBarrier * aCoulombBarrier) 49 : theA(anA), 50 theZ(aZ), 51 Gamma(aGamma), 52 theCoulombBarrierptr(aCoulombBarrier) 53 {} 54 55 G4EvaporationProbability::G4EvaporationProbability() 56 : theA(0), 57 theZ(0), 58 Gamma(0.0), 59 theCoulombBarrierptr(0) 60 {} 61 62 G4EvaporationProbability::~G4EvaporationProbability() 63 {} 64 65 G4double 66 G4EvaporationProbability::EmissionProbability(const G4Fragment & fragment, G4double anEnergy) 48 67 { 49 throw G4HadronicException(__FILE__, __LINE__, "G4EvaporationProbability::copy_constructor meant to not be accessable"); 50 } 51 52 53 54 55 const G4EvaporationProbability & G4EvaporationProbability:: 56 operator=(const G4EvaporationProbability &) 68 G4double EmissionProbability = 0.0; 69 G4double MaximalKineticEnergy = anEnergy; 70 71 if (MaximalKineticEnergy > 0.0 && fragment.GetExcitationEnergy() > 0.0) { 72 EmissionProbability = CalculateProbability(fragment, MaximalKineticEnergy); 73 74 } 75 return EmissionProbability; 76 } 77 78 //////////////////////////////////// 79 80 // Computes the integrated probability of evaporation channel 81 G4double 82 G4EvaporationProbability::CalculateProbability(const G4Fragment & fragment, 83 G4double MaximalKineticEnergy) 57 84 { 58 throw G4HadronicException(__FILE__, __LINE__, "G4EvaporationProbability::operator= meant to not be accessable"); 59 return *this; 60 } 61 62 63 G4bool G4EvaporationProbability::operator==(const G4EvaporationProbability &) const 64 { 65 return false; 66 } 67 68 G4bool G4EvaporationProbability::operator!=(const G4EvaporationProbability &) const 69 { 70 return true; 71 } 72 73 G4double G4EvaporationProbability::EmissionProbability(const G4Fragment & fragment, const G4double anEnergy) 74 { 75 G4double EmissionProbability = 0.0; 76 G4double MaximalKineticEnergy = anEnergy; 77 78 if (MaximalKineticEnergy > 0.0 && fragment.GetExcitationEnergy() > 0.0) { 79 EmissionProbability = CalculateProbability(fragment, MaximalKineticEnergy); 80 81 } 82 return EmissionProbability; 83 } 84 85 //////////////////////////////////// 86 87 // Computes the integrated probability of evaporation channel 88 G4double G4EvaporationProbability::CalculateProbability(const G4Fragment & fragment, const G4double MaximalKineticEnergy) 89 { 90 G4double ResidualA = fragment.GetA() - theA; 91 G4double ResidualZ = fragment.GetZ() - theZ; 92 G4double U = fragment.GetExcitationEnergy(); 85 G4int ResidualA = fragment.GetA_asInt() - theA; 86 G4int ResidualZ = fragment.GetZ_asInt() - theZ; 87 G4double U = fragment.GetExcitationEnergy(); 93 88 94 if (OPTxs==0) { 95 96 97 G4double NuclearMass = G4ParticleTable::GetParticleTable()->GetIonTable()->GetNucleusMass(theZ,theA); 98 99 100 G4double delta0 = G4PairingCorrection::GetInstance()->GetPairingCorrection(static_cast<G4int>(fragment.GetA()), 101 static_cast<G4int>(fragment.GetZ())); 102 103 G4double SystemEntropy = 2.0*std::sqrt(theEvapLDPptr->LevelDensityParameter(static_cast<G4int>(fragment.GetA()), 104 static_cast<G4int>(fragment.GetZ()),U)* 105 (U-delta0)); 89 if (OPTxs==0) { 90 91 G4double NuclearMass = fragment.ComputeGroundStateMass(theZ,theA); 92 93 G4double delta0 = fPairCorr->GetPairingCorrection(fragment.GetA_asInt(), 94 fragment.GetZ_asInt()); 95 96 G4double SystemEntropy = 2.0*std::sqrt( 97 theEvapLDPptr->LevelDensityParameter(fragment.GetA_asInt(),fragment.GetZ_asInt(),U)* 98 (U-delta0)); 106 99 107 108 G4double RN = 1.5*fermi; 100 const G4double RN = 1.5*fermi; 109 101 110 102 G4double Alpha = CalcAlphaParam(fragment); … … 112 104 113 105 G4double Rmax = MaximalKineticEnergy; 114 G4double a = theEvapLDPptr->LevelDensityParameter(static_cast<G4int>(ResidualA), 115 static_cast<G4int>(ResidualZ), 116 Rmax); 117 G4double GlobalFactor = static_cast<G4double>(Gamma) * (Alpha/(a*a)) * 118 (NuclearMass*RN*RN*std::pow(ResidualA,2./3.))/ 119 (2.*pi* hbar_Planck*hbar_Planck); 106 G4double a = theEvapLDPptr->LevelDensityParameter(ResidualA,ResidualZ,Rmax); 107 G4double GlobalFactor = Gamma * Alpha/(a*a) * 108 (NuclearMass*RN*RN*fG4pow->Z23(ResidualA))/ 109 (twopi* hbar_Planck*hbar_Planck); 120 110 G4double Term1 = (2.0*Beta*a-3.0)/2.0 + Rmax*a; 121 111 G4double Term2 = (2.0*Beta*a-3.0)*std::sqrt(Rmax*a) + 2.0*a*Rmax; 122 112 123 113 G4double ExpTerm1 = 0.0; 124 if (SystemEntropy <= 600.0) ExpTerm1 = std::exp(-SystemEntropy);114 if (SystemEntropy <= 600.0) { ExpTerm1 = std::exp(-SystemEntropy); } 125 115 126 116 G4double ExpTerm2 = 2.*std::sqrt(a*Rmax) - SystemEntropy; 127 if (ExpTerm2 > 700.0) ExpTerm2 = 700.0;117 if (ExpTerm2 > 700.0) { ExpTerm2 = 700.0; } 128 118 ExpTerm2 = std::exp(ExpTerm2); 129 119 … … 134 124 } else if (OPTxs==1 || OPTxs==2 ||OPTxs==3 || OPTxs==4) { 135 125 136 G4double limit; 137 if (useSICB) limit=theCoulombBarrierptr->GetCoulombBarrier(G4lrint(ResidualA),G4lrint(ResidualZ),U); 138 else limit=0.; 139 140 if (MaximalKineticEnergy <= limit) return 0.0; 141 142 143 // if Coulomb barrier cutoff is superimposed for all cross sections the limit is the Coulomb Barrier 126 G4double EvaporatedMass = fragment.ComputeGroundStateMass(theZ,theA); 127 G4double ResidulalMass = fragment.ComputeGroundStateMass(ResidualZ,ResidualA); 128 G4double limit = std::max(0.0,fragment.GetGroundStateMass()-EvaporatedMass-ResidulalMass); 129 if (useSICB) { 130 limit = std::max(limit,theCoulombBarrierptr->GetCoulombBarrier(ResidualA,ResidualZ,U)); 131 } 132 133 if (MaximalKineticEnergy <= limit) { return 0.0; } 134 135 // if Coulomb barrier cutoff is superimposed for all cross sections 136 // then the limit is the Coulomb Barrier 144 137 G4double LowerLimit= limit; 145 138 146 // 139 //MaximalKineticEnergy: asimptotic value (already accounted for in G4EvaporationChannel) 147 140 148 141 G4double UpperLimit = MaximalKineticEnergy; 149 142 150 151 143 G4double Width = IntegrateEmissionProbability(fragment,LowerLimit,UpperLimit); 152 144 153 145 return Width; 154 } else {146 } else { 155 147 std::ostringstream errOs; 156 148 errOs << "Bad option for cross sections at evaporation" <<G4endl; … … 163 155 164 156 G4double G4EvaporationProbability:: 165 IntegrateEmissionProbability(const G4Fragment & fragment, const G4double & Low, const G4double & Up ) 157 IntegrateEmissionProbability(const G4Fragment & fragment, 158 const G4double & Low, const G4double & Up ) 166 159 { 167 168 160 static const G4int N = 10; 169 161 // 10-Points Gauss-Legendre abcisas and weights … … 212 204 //New method (OPT=1,2,3,4) 213 205 214 G4double G4EvaporationProbability::ProbabilityDistributionFunction( const G4Fragment & fragment, const G4double K) 206 G4double 207 G4EvaporationProbability::ProbabilityDistributionFunction( const G4Fragment & fragment, 208 G4double K) 215 209 { 216 210 G4int ResidualA = fragment.GetA_asInt() - theA; 211 G4int ResidualZ = fragment.GetZ_asInt() - theZ; 212 G4double U = fragment.GetExcitationEnergy(); 213 //G4cout << "### G4EvaporationProbability::ProbabilityDistributionFunction" << G4endl; 214 //G4cout << "FragZ= " << fragment.GetZ_asInt() << " FragA= " << fragment.GetA_asInt() 215 // << " Z= " << theZ << " A= " << theA << G4endl; 216 //G4cout << "PC " << fPairCorr << " DP " << theEvapLDPptr << G4endl; 217 218 // if(K <= theCoulombBarrierptr->GetCoulombBarrier(ResidualA,ResidualZ,U)) return 0.0; 219 220 G4double delta1 = fPairCorr->GetPairingCorrection(ResidualA,ResidualZ); 217 221 218 219 220 G4double ResidualA = fragment.GetA() - theA; 221 G4double ResidualZ = fragment.GetZ() - theZ; 222 G4double U = fragment.GetExcitationEnergy(); 223 224 // if(K <= theCoulombBarrierptr->GetCoulombBarrier(G4lrint(ResidualA),G4lrint(ResidualZ),U)) return 0.0; 225 226 G4double delta1 = G4PairingCorrection::GetInstance()->GetPairingCorrection(static_cast<G4int>(ResidualA),static_cast<G4int>(ResidualZ)); 227 228 229 G4double delta0 = G4PairingCorrection::GetInstance()->GetPairingCorrection(static_cast<G4int>(fragment.GetA()),static_cast<G4int>(fragment.GetZ())); 230 231 232 G4double ParticleMass = G4ParticleTable::GetParticleTable()->GetIonTable()->GetNucleusMass(theZ,theA); 233 234 G4double theSeparationEnergy= G4NucleiProperties::GetMassExcess(static_cast<G4int>(theA),static_cast<G4int>(theZ)) + 235 G4NucleiProperties::GetMassExcess(static_cast<G4int>(ResidualA),static_cast<G4int>(ResidualZ)) - 236 G4NucleiProperties::GetMassExcess(static_cast<G4int>(fragment.GetA()),static_cast<G4int>(fragment.GetZ())); 237 238 G4double a0 = theEvapLDPptr->LevelDensityParameter(static_cast<G4int>(fragment.GetA()),static_cast<G4int>(fragment.GetZ()),U - delta0); 239 240 G4double a1 = theEvapLDPptr->LevelDensityParameter(static_cast<G4int>(ResidualA),static_cast<G4int>(ResidualZ),U - theSeparationEnergy - delta1); 241 242 243 G4double E0=U-delta0; 244 245 G4double E1=U-theSeparationEnergy-delta1-K; 246 247 if (E1<0.) return 0.; 222 G4double delta0 = fPairCorr->GetPairingCorrection(fragment.GetA_asInt(), 223 fragment.GetZ_asInt()); 224 225 226 G4double ParticleMass = fragment.ComputeGroundStateMass(theZ,theA); 227 G4double ResidualMass = fragment.ComputeGroundStateMass(ResidualZ,ResidualA); 228 229 G4double theSeparationEnergy = ParticleMass + ResidualMass 230 - fragment.GetGroundStateMass(); 231 232 G4double a0 = theEvapLDPptr->LevelDensityParameter(fragment.GetA_asInt(), 233 fragment.GetZ_asInt(), 234 U - delta0); 235 236 G4double a1 = theEvapLDPptr->LevelDensityParameter(ResidualA, ResidualZ, 237 U - theSeparationEnergy - delta1); 238 239 240 G4double E0 = U - delta0; 241 242 G4double E1 = U - theSeparationEnergy - delta1 - K; 243 244 if (E1<0.) { return 0.; } 248 245 249 246 //JMQ 14/02/09 BUG fixed: hbarc should be in the denominator instead of hbar_Planck 250 247 //Without 1/hbar_Panck remains as a width 251 // G4double Prob=Gamma*ParticleMass/((pi*hbar_Planck)*(pi*hbar_Planck)*252 //std::exp(2*std::sqrt(a0*E0)))*K*CrossSection(fragment,K)*std::exp(2*std::sqrt(a1*E1))*millibarn;253 248 254 249 G4double Prob=Gamma*ParticleMass/((pi*hbarc)*(pi*hbarc)*std::exp(2*std::sqrt(a0*E0))) -
trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4He3EvaporationChannel.cc
r1340 r1347 25 25 // 26 26 // 27 // $Id: G4He3EvaporationChannel.cc,v 1. 4 2006/06/29 20:10:33 gunterExp $28 // GEANT4 tag $Name: geant4-09-0 3-ref-09$27 // $Id: G4He3EvaporationChannel.cc,v 1.5 2010/11/17 12:14:59 vnivanch Exp $ 28 // GEANT4 tag $Name: geant4-09-04-ref-00 $ 29 29 // 30 30 // Hadronic Process: Nuclear De-excitations 31 31 // by V. Lara (Nov. 1999) 32 32 // 33 // 17-11-2010 V.Ivanchenko moved constructor and destructor to source and cleanup 33 34 34 35 #include "G4He3EvaporationChannel.hh" 35 36 37 G4He3EvaporationChannel::G4He3EvaporationChannel() 38 : G4EvaporationChannel(3,2,"He3",&theEvaporationProbability,&theCoulombBarrier) 39 {} 36 40 37 const G4He3EvaporationChannel & G4He3EvaporationChannel::operator=(const G4He3EvaporationChannel & ) 38 { 39 throw G4HadronicException(__FILE__, __LINE__, "G4He3EvaporationChannel::operator= meant to not be accessable"); 40 return *this; 41 } 41 G4He3EvaporationChannel::~G4He3EvaporationChannel() 42 {} 42 43 43 G4He3EvaporationChannel::G4He3EvaporationChannel(const G4He3EvaporationChannel & ) : G4EvaporationChannel()44 {45 throw G4HadronicException(__FILE__, __LINE__, "G4He3EvaporationChannel::CopyConstructor meant to not be accessable");46 }47 48 G4bool G4He3EvaporationChannel::operator==(const G4He3EvaporationChannel & right) const49 {50 return (this == (G4He3EvaporationChannel *) &right);51 // return false;52 }53 54 G4bool G4He3EvaporationChannel::operator!=(const G4He3EvaporationChannel & right) const55 {56 return (this != (G4He3EvaporationChannel *) &right);57 // return true;58 }59 -
trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4He3EvaporationProbability.cc
r1315 r1347 24 24 // ******************************************************************** 25 25 // 26 //J.M. Quesada (August2008). Based on: 26 // $Id: G4He3EvaporationProbability.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 "G4He3EvaporationProbability.hh" … … 36 40 G4He3EvaporationProbability::G4He3EvaporationProbability() : 37 41 G4EvaporationProbability(3,2,2,&theCoulombBarrier) // A,Z,Gamma,&theCoulombBarrier 38 { 39 40 } 41 42 G4He3EvaporationProbability::G4He3EvaporationProbability(const G4He3EvaporationProbability &) : G4EvaporationProbability() 43 { 44 throw G4HadronicException(__FILE__, __LINE__, "G4He3EvaporationProbability::copy_constructor meant to not be accessable"); 45 } 46 47 48 49 50 const G4He3EvaporationProbability & G4He3EvaporationProbability:: 51 operator=(const G4He3EvaporationProbability &) 52 { 53 throw G4HadronicException(__FILE__, __LINE__, "G4He3EvaporationProbability::operator= meant to not be accessable"); 54 return *this; 55 } 56 57 58 G4bool G4He3EvaporationProbability::operator==(const G4He3EvaporationProbability &) const 59 { 60 return false; 61 } 62 63 G4bool G4He3EvaporationProbability::operator!=(const G4He3EvaporationProbability &) const 64 { 65 return true; 66 } 67 68 G4double G4He3EvaporationProbability::CalcAlphaParam(const G4Fragment & fragment) 69 { return 1.0 + CCoeficient(static_cast<G4double>(fragment.GetZ()-GetZ()));} 42 {} 43 44 G4He3EvaporationProbability::~G4He3EvaporationProbability() 45 {} 46 47 G4double G4He3EvaporationProbability::CalcAlphaParam(const G4Fragment & fragment) 48 { return 1.0 + CCoeficient(fragment.GetZ_asInt()-GetZ());} 70 49 71 50 G4double G4He3EvaporationProbability::CalcBetaParam(const G4Fragment & ) 72 51 { return 0.0; } 73 52 74 53 75 G4double G4He3EvaporationProbability::CCoeficient( const G4doubleaZ)76 { 77 78 79 80 81 82 83 84 85 54 G4double G4He3EvaporationProbability::CCoeficient(G4int aZ) 55 { 56 // Data comes from 57 // Dostrovsky, Fraenkel and Friedlander 58 // Physical Review, vol 116, num. 3 1959 59 // 60 // const G4int size = 5; 61 // G4double Zlist[5] = { 10.0, 20.0, 30.0, 50.0, 70.0}; 62 // G4double Calpha[5] = { 0.10, 0.10, 0.10, 0.08, 0.06}; 63 // C for He3 is equal to C for alpha times 4/3 64 G4double C = 0.0; 86 65 87 88 if (aZ <= 30) { 89 C = 0.10; 90 } else if (aZ <= 50) { 91 C = 0.1 + -((aZ-50.)/20.)*0.02; 92 } else if (aZ < 70) { 93 C = 0.08 + -((aZ-70.)/20.)*0.02; 94 } else { 95 C = 0.06; 96 } 97 return C*(4.0/3.0); 66 if (aZ <= 30) 67 { 68 C = 0.10; 69 } 70 else if (aZ <= 50) 71 { 72 C = 0.1 - (aZ - 30)*0.001; 73 } 74 else if (aZ < 70) 75 { 76 C = 0.08 - (aZ - 50)*0.001; 77 } 78 else 79 { 80 C = 0.06; 81 } 82 return C*(4.0/3.0); 98 83 } 99 84 … … 104 89 //OPT=3,4 Kalbach's parameterization 105 90 // 106 G4double 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 ************************ 91 G4double 92 G4He3EvaporationProbability::CrossSection(const G4Fragment & fragment, G4double K) 93 { 94 95 theA=GetA(); 96 theZ=GetZ(); 97 ResidualA=fragment.GetA_asInt()-theA; 98 ResidualZ=fragment.GetZ_asInt()-theZ; 99 100 ResidualAthrd=fG4pow->Z13(ResidualA); 101 FragmentA=fragment.GetA_asInt(); 102 FragmentAthrd=fG4pow->Z13(FragmentA); 103 104 if (OPTxs==0) {std::ostringstream errOs; 105 errOs << "We should'n be here (OPT =0) at evaporation cross section calculation (He3's)!!" 106 <<G4endl; 107 throw G4HadronicException(__FILE__, __LINE__, errOs.str()); 108 return 0.;} 109 if( OPTxs==1 || OPTxs==2) return G4He3EvaporationProbability::GetOpt12( K); 110 else if (OPTxs==3 || OPTxs==4) return G4He3EvaporationProbability::GetOpt34( K); 111 else{ 112 std::ostringstream errOs; 113 errOs << "BAD He3's CROSS SECTION OPTION AT EVAPORATION!!" <<G4endl; 114 throw G4HadronicException(__FILE__, __LINE__, errOs.str()); 115 return 0.; 116 } 117 } 118 119 //********************* OPT=1,2 : Chatterjee's cross section ***************** 134 120 //(fitting to cross section from Bechetti & Greenles OM potential) 135 121 136 122 G4double G4He3EvaporationProbability::GetOpt12(const G4double K) 137 123 { 138 139 G4double 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;} 124 G4double Kc = K; 125 126 // JMQ xsec is set constat above limit of validity 127 if (K > 50*MeV) { Kc = 50*MeV; } 128 129 G4double landa ,mu ,nu ,p , Ec,q,r,ji,xs; 130 131 G4double p0 = -3.06; 132 G4double p1 = 278.5; 133 G4double p2 = -1389.; 134 G4double landa0 = -0.00535; 135 G4double landa1 = -11.16; 136 G4double mu0 = 555.5; 137 G4double mu1 = 0.40; 138 G4double nu0 = 687.4; 139 G4double nu1 = -476.3; 140 G4double nu2 = 0.509; 141 G4double delta=1.2; 142 143 Ec = 1.44*theZ*ResidualZ/(1.5*ResidualAthrd+delta); 144 p = p0 + p1/Ec + p2/(Ec*Ec); 145 landa = landa0*ResidualA + landa1; 146 147 G4double resmu1 = fG4pow->powZ(ResidualA,mu1); 148 mu = mu0*resmu1; 149 nu = resmu1*(nu0 + nu1*Ec + nu2*(Ec*Ec)); 150 q = landa - nu/(Ec*Ec) - 2*p*Ec; 151 r = mu + 2*nu/Ec + p*(Ec*Ec); 152 153 ji=std::max(Kc,Ec); 154 if(Kc < Ec) { xs = p*Kc*Kc + q*Kc + r;} 155 else {xs = p*(Kc - ji)*(Kc - ji) + landa*Kc + mu + nu*(2 - Kc/ji)/ji ;} 156 157 if (xs <0.0) {xs=0.0;} 171 158 172 159 return xs; 173 160 174 161 } … … 178 165 //c ** 3he from o.m. of gibson et al 179 166 { 180 181 167 G4double landa, mu, nu, p , signor(1.),sig; 182 168 G4double ec,ecsq,xnulam,etest(0.),a; 183 169 G4double b,ecut,cut,ecut2,geom,elab; 184 185 170 186 171 G4double flow = 1.e-18; 187 172 G4double spill= 1.e+18; 188 189 173 190 174 G4double p0 = -2.88; … … 200 184 201 185 G4double ra=0.80; 202 186 203 187 //JMQ 13/02/09 increase of reduced radius to lower the barrier 204 188 // ec = 1.44 * theZ * ResidualZ / (1.5*ResidualAthrd+ra); … … 207 191 p = p0 + p1/ec + p2/ecsq; 208 192 landa = landa0*ResidualA + landa1; 209 a = std::pow(ResidualA,mu1);193 a = fG4pow->powZ(ResidualA,mu1); 210 194 mu = mu0 * a; 211 195 nu = a* (nu0+nu1*ec+nu2*ecsq); 212 196 xnulam = nu / landa; 213 if (xnulam > spill) xnulam=0.;214 if (xnulam >= flow) etest = 1.2 *std::sqrt(xnulam);197 if (xnulam > spill) { xnulam=0.; } 198 if (xnulam >= flow) { etest = 1.2 *std::sqrt(xnulam); } 215 199 216 200 a = -2.*p*ec + landa - nu/ecsq; … … 221 205 ecut = (ecut-a) / (p+p); 222 206 ecut2 = ecut; 223 //JMQ 290310 for avoiding unphysical increase below minimum (at ecut) 224 //ecut<0 means that there is no cut with energy axis, i.e. xs is set to 0 bellow minimum 225 // if (cut < 0.) ecut2 = ecut - 2.; 226 if (cut < 0.) ecut2 = ecut; 227 elab = K * FragmentA / ResidualA; 207 //JMQ 290310 for avoiding unphysical increase below minimum (at ecut) 208 // ecut<0 means that there is no cut with energy axis, i.e. xs is set 209 // to 0 bellow minimum 210 // if (cut < 0.) ecut2 = ecut - 2.; 211 if (cut < 0.) { ecut2 = ecut; } 212 elab = K * FragmentA /G4double(ResidualA); 228 213 sig = 0.; 229 214 230 215 if (elab <= ec) { //start for E<Ec 231 if (elab > ecut2) sig = (p*elab*elab+a*elab+b) * signor;216 if (elab > ecut2) { sig = (p*elab*elab+a*elab+b) * signor; } 232 217 } //end for E<Ec 233 218 else { //start for E>Ec 234 219 sig = (landa*elab+mu+nu/elab) * signor; 235 220 geom = 0.; 236 if (xnulam < flow || elab < etest) return sig;221 if (xnulam < flow || elab < etest) { return sig; } 237 222 geom = std::sqrt(theA*K); 238 223 geom = 1.23*ResidualAthrd + ra + 4.573/geom; … … 243 228 244 229 } 245 246 // ************************** end of cross sections *******************************247 -
trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4NeutronEvaporationChannel.cc
r1340 r1347 25 25 // 26 26 // 27 // $Id: G4NeutronEvaporationChannel.cc,v 1. 4 2006/06/29 20:10:37 gunterExp $28 // GEANT4 tag $Name: geant4-09-0 3-ref-09$27 // $Id: G4NeutronEvaporationChannel.cc,v 1.5 2010/11/17 12:14:59 vnivanch Exp $ 28 // GEANT4 tag $Name: geant4-09-04-ref-00 $ 29 29 // 30 30 // Hadronic Process: Nuclear De-excitations 31 31 // by V. Lara (Nov. 1999) 32 32 // 33 // 17-11-2010 V.Ivanchenko moved constructor and destructor to source and cleanup 33 34 34 35 #include "G4NeutronEvaporationChannel.hh" 35 36 37 G4NeutronEvaporationChannel::G4NeutronEvaporationChannel() 38 : G4EvaporationChannel(1,0,"neutron",&theEvaporationProbability,&theCoulombBarrier) 39 {} 36 40 37 const G4NeutronEvaporationChannel & G4NeutronEvaporationChannel::operator=(const G4NeutronEvaporationChannel & ) 38 { 39 throw G4HadronicException(__FILE__, __LINE__, "G4NeutronEvaporationChannel::operator= meant to not be accessable"); 40 return *this; 41 } 42 43 G4NeutronEvaporationChannel::G4NeutronEvaporationChannel(const G4NeutronEvaporationChannel & ) : G4EvaporationChannel() 44 { 45 throw G4HadronicException(__FILE__, __LINE__, "G4NeutronEvaporationChannel::CopyConstructor meant to not be accessable"); 46 } 47 48 G4bool G4NeutronEvaporationChannel::operator==(const G4NeutronEvaporationChannel & right) const 49 { 50 return (this == (G4NeutronEvaporationChannel *) &right); 51 // return false; 52 } 53 54 G4bool G4NeutronEvaporationChannel::operator!=(const G4NeutronEvaporationChannel & right) const 55 { 56 return (this != (G4NeutronEvaporationChannel *) &right); 57 // return true; 58 } 59 41 G4NeutronEvaporationChannel::~G4NeutronEvaporationChannel() 42 {} -
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 -
trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4ProtonEvaporationChannel.cc
r1340 r1347 25 25 // 26 26 // 27 // $Id: G4ProtonEvaporationChannel.cc,v 1. 4 2006/06/29 20:10:41 gunterExp $28 // GEANT4 tag $Name: geant4-09-0 3-ref-09$27 // $Id: G4ProtonEvaporationChannel.cc,v 1.5 2010/11/17 12:14:59 vnivanch Exp $ 28 // GEANT4 tag $Name: geant4-09-04-ref-00 $ 29 29 // 30 30 // Hadronic Process: Nuclear De-excitations 31 31 // by V. Lara (Nov. 1999) 32 32 // 33 // 17-11-2010 V.Ivanchenko moved constructor and destructor to source and cleanup 33 34 34 35 #include "G4ProtonEvaporationChannel.hh" 35 36 37 G4ProtonEvaporationChannel::G4ProtonEvaporationChannel() 38 : G4EvaporationChannel(1,1,"proton",&theEvaporationProbability,&theCoulombBarrier) 39 {} 36 40 37 const G4ProtonEvaporationChannel & G4ProtonEvaporationChannel::operator=(const G4ProtonEvaporationChannel & ) 38 { 39 throw G4HadronicException(__FILE__, __LINE__, "G4ProtonEvaporationChannel::operator= meant to not be accessable"); 40 return *this; 41 } 42 43 G4ProtonEvaporationChannel::G4ProtonEvaporationChannel(const G4ProtonEvaporationChannel & ) : G4EvaporationChannel() 44 { 45 throw G4HadronicException(__FILE__, __LINE__, "G4ProtonEvaporationChannel::CopyConstructor meant to not be accessable"); 46 } 47 48 G4bool G4ProtonEvaporationChannel::operator==(const G4ProtonEvaporationChannel & right) const 49 { 50 return (this == (G4ProtonEvaporationChannel *) &right); 51 // return false; 52 } 53 54 G4bool G4ProtonEvaporationChannel::operator!=(const G4ProtonEvaporationChannel & right) const 55 { 56 return (this != (G4ProtonEvaporationChannel *) &right); 57 // return true; 58 } 59 60 41 G4ProtonEvaporationChannel::~G4ProtonEvaporationChannel() 42 {} -
trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4ProtonEvaporationProbability.cc
r1315 r1347 24 24 // ******************************************************************** 25 25 // 26 //J.M. Quesada (August2008). Based on: 26 // $Id: G4ProtonEvaporationProbability.cc,v 1.17 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 "G4ProtonEvaporationProbability.hh" … … 36 40 G4ProtonEvaporationProbability::G4ProtonEvaporationProbability() : 37 41 G4EvaporationProbability(1,1,2,&theCoulombBarrier) // A,Z,Gamma,&theCoulombBarrier 38 { 39 40 } 41 42 G4ProtonEvaporationProbability::G4ProtonEvaporationProbability(const G4ProtonEvaporationProbability &) : G4EvaporationProbability() 43 { 44 throw G4HadronicException(__FILE__, __LINE__, "G4ProtonEvaporationProbability::copy_constructor meant to not be accessable"); 45 } 46 47 const G4ProtonEvaporationProbability & G4ProtonEvaporationProbability:: 48 operator=(const G4ProtonEvaporationProbability &) 49 { 50 throw G4HadronicException(__FILE__, __LINE__, "G4ProtonEvaporationProbability::operator= meant to not be accessable"); 51 return *this; 52 } 53 54 55 G4bool G4ProtonEvaporationProbability::operator==(const G4ProtonEvaporationProbability &) const 56 { 57 return false; 58 } 59 60 G4bool G4ProtonEvaporationProbability::operator!=(const G4ProtonEvaporationProbability &) const 61 { 62 return true; 63 } 64 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 & ) 42 {} 43 44 G4ProtonEvaporationProbability::~G4ProtonEvaporationProbability() 45 {} 46 47 G4double G4ProtonEvaporationProbability::CalcAlphaParam(const G4Fragment & fragment) 48 { return 1.0 + CCoeficient(fragment.GetZ_asInt()-GetZ());} 49 50 G4double G4ProtonEvaporationProbability::CalcBetaParam(const G4Fragment & ) 69 51 { return 0.0; } 70 52 71 G4double G4ProtonEvaporationProbability::CCoeficient(const G4doubleaZ)72 { 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 53 G4double G4ProtonEvaporationProbability::CCoeficient(G4int aZ) 54 { 55 // Data comes from 56 // Dostrovsky, Fraenkel and Friedlander 57 // Physical Review, vol 116, num. 3 1959 58 // 59 // const G4int size = 5; 60 // G4double Zlist[5] = { 10.0, 20.0, 30.0, 50.0, 70.0}; 61 // G4double Cp[5] = { 0.50, 0.28, 0.20, 0.15, 0.10}; 62 G4double C = 0.0; 63 64 if (aZ >= 70) { 65 C = 0.10; 66 } else { 67 C = ((((0.15417e-06*aZ) - 0.29875e-04)*aZ + 0.21071e-02)*aZ - 0.66612e-01)*aZ + 0.98375; 68 } 69 70 return C; 89 71 90 72 } … … 97 79 //OPT=3 Kalbach's parameterization 98 80 // 99 G4double G4ProtonEvaporationProbability::CrossSection(const G4Fragment & fragment, const G4double K) 100 { 101 // G4cout<<" In G4ProtonEVaporationProbability OPTxs="<<OPTxs<<G4endl; 102 // G4cout<<" In G4ProtonEVaporationProbability useSICB="<<useSICB<<G4endl; 81 G4double 82 G4ProtonEvaporationProbability::CrossSection(const G4Fragment & fragment, G4double K) 83 { 84 // G4cout<<" In G4ProtonEVaporationProbability OPTxs="<<OPTxs<<G4endl; 85 // G4cout<<" In G4ProtonEVaporationProbability useSICB="<<useSICB<<G4endl; 103 86 104 87 theA=GetA(); 105 88 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); 89 ResidualA=fragment.GetA_asInt()-theA; 90 ResidualZ=fragment.GetZ_asInt()-theZ; 91 92 ResidualAthrd=fG4pow->Z13(ResidualA); 93 FragmentA=fragment.GetA_asInt(); 94 FragmentAthrd=fG4pow->Z13(FragmentA); 95 112 96 U=fragment.GetExcitationEnergy(); 113 97 … … 126 110 } 127 111 } 128 //********************* OPT=1 : Chatterjee's cross section ************************ 112 113 //********************* OPT=1 : Chatterjee's cross section ********************* 129 114 //(fitting to cross section from Bechetti & Greenles OM potential) 130 115 131 G4double G4ProtonEvaporationProbability::GetOpt1( constG4double K)116 G4double G4ProtonEvaporationProbability::GetOpt1(G4double K) 132 117 { 133 118 G4double Kc=K; 134 119 135 // JMQ xsec is set constat above limit of validity136 if (K >50) Kc=50;120 // JMQ xsec is set constat above limit of validity 121 if (K > 50*MeV) { Kc = 50*MeV; } 137 122 138 123 G4double landa, landa0, landa1, mu, mu0, mu1,nu, nu0, nu1, nu2,xs; … … 154 139 p = p0 + p1/Ec + p2/(Ec*Ec); 155 140 landa = landa0*ResidualA + landa1; 156 mu = mu0*std::pow(ResidualA,mu1); 157 nu = std::pow(ResidualA,mu1)*(nu0 + nu1*Ec + nu2*(Ec*Ec)); 141 142 G4double resmu1 = fG4pow->powZ(ResidualA,mu1); 143 mu = mu0*resmu1; 144 nu = resmu1*(nu0 + nu1*Ec + nu2*(Ec*Ec)); 158 145 q = landa - nu/(Ec*Ec) - 2*p*Ec; 159 146 r = mu + 2*nu/Ec + p*(Ec*Ec); … … 162 149 if(Kc < Ec) { xs = p*Kc*Kc + q*Kc + r;} 163 150 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 174 G4double G4ProtonEvaporationProbability::GetOpt2(const G4double K) 175 { 176 177 G4double rnpro,rnneu,eekin,ekin,ff1,ff2,ff3,r0,fac,fac1,fac2,b0,xine_th(0); 151 if (xs <0.0) {xs=0.0;} 152 153 return xs; 154 } 155 156 //************* OPT=2 : Welisch's proton reaction cross section *************** 157 158 G4double G4ProtonEvaporationProbability::GetOpt2(G4double K) 159 { 160 161 G4double eekin,ekin,ff1,ff2,ff3,r0,fac,fac1,fac2,b0,xine_th(0); 178 162 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; 163 // This is redundant when the Coulomb barrier is overimposed to all 164 // cross sections 165 // It should be kept when Coulomb barrier only imposed at OPTxs=2 166 167 if(!useSICB && K<=theCoulombBarrier.GetCoulombBarrier(ResidualA,ResidualZ,U)) 168 { return 0.0; } 183 169 184 170 eekin=K; 185 rnpro=ResidualZ; 186 rnneu=ResidualA-ResidualZ; 171 G4int rnneu=ResidualA-ResidualZ; 187 172 ekin=eekin/1000; 188 189 173 r0=1.36*1.e-15; 190 174 fac=pi*r0*r0; … … 192 176 fac1=b0*(1.-1./ResidualAthrd); 193 177 fac2=1.; 194 if(rnneu > 1.5) fac2=std::log(rnneu);178 if(rnneu > 1.5) { fac2 = fG4pow->logZ(rnneu); } 195 179 xine_th= 1.e+31*fac*fac2*(1.+ResidualAthrd-fac1); 196 180 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;181 ff1=0.70-0.0020*ResidualA; 182 ff2=1.00+1/G4double(ResidualA); 183 ff3=0.8+18/G4double(ResidualA)-0.002*ResidualA; 200 184 fac=1.-(1./(1.+std::exp(-8.*ff1*(std::log10(ekin)+1.37*ff2)))); 201 185 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;186 ff1=1.-1/G4double(ResidualA)-0.001*ResidualA; 187 ff2=1.17-2.7/G4double(ResidualA)-0.0014*ResidualA; 204 188 fac=-8.*ff1*(std::log10(ekin)+2.0*ff2); 205 189 fac=1./(1.+std::exp(fac)); 206 xine_th=xine_th*fac; 190 xine_th=xine_th*fac; 207 191 if (xine_th < 0.0){ 208 192 std::ostringstream errOs; … … 212 196 throw G4HadronicException(__FILE__, __LINE__, errOs.str()); 213 197 } 214 215 198 return xine_th; 216 217 } 218 199 } 219 200 220 201 // *********** OPT=3 : Kalbach's cross sections (from PRECO code)************* 221 202 G4double G4ProtonEvaporationProbability::GetOpt3(const G4double K) 222 203 { 223 // ** p from becchetti and greenlees (but modified with sub-barrier224 // ** correction function and xp2 changed from -449)204 // ** p from becchetti and greenlees (but modified with sub-barrier 205 // ** correction function and xp2 changed from -449) 225 206 226 207 G4double landa, landa0, landa1, mu, mu0, mu1,nu, nu0, nu1, nu2; … … 236 217 nu1 = -182.4; 237 218 nu2 = -1.872; 238 239 // parameters for proton cross section refinement219 220 // parameters for proton cross section refinement 240 221 G4double afit,bfit,a2,b2; 241 222 afit=-0.0785656; … … 243 224 a2= -0.00089076; 244 225 b2= 0.0231597; 245 226 246 227 G4double ec,ecsq,xnulam,etest(0.),ra(0.),a,w,c,signor(1.),signor2,sig; 247 228 G4double b,ecut,cut,ecut2,geom,elab; 248 249 229 250 230 G4double flow = 1.e-18; 251 231 G4double spill= 1.e+18; 252 253 254 if (ResidualA <= 60.) signor = 0.92; 255 else if (ResidualA < 100.) signor = 0.8 + ResidualA*0.002; 256 257 232 233 if (ResidualA <= 60) { signor = 0.92; } 234 else if (ResidualA < 100) { signor = 0.8 + ResidualA*0.002; } 235 258 236 ec = 1.44 * theZ * ResidualZ / (1.5*ResidualAthrd+ra); 259 237 ecsq = ec * ec; 260 238 p = p0 + p1/ec + p2/ecsq; 261 239 landa = landa0*ResidualA + landa1; 262 a = std::pow(ResidualA,mu1);240 a = fG4pow->powZ(ResidualA,mu1); 263 241 mu = mu0 * a; 264 242 nu = a* (nu0+nu1*ec+nu2*ecsq); 265 243 266 244 c =std::min(3.15,ec*0.5); 267 245 w = 0.7 * c / 3.15; 268 246 269 247 xnulam = nu / landa; 270 if (xnulam > spill) xnulam=0.;271 if (xnulam >= flow) etest =std::sqrt(xnulam) + 7.;272 248 if (xnulam > spill) { xnulam=0.; } 249 if (xnulam >= flow) { etest =std::sqrt(xnulam) + 7.; } 250 273 251 a = -2.*p*ec + landa - nu/ecsq; 274 252 b = p*ecsq + mu + 2.*nu/ec; 275 253 ecut = 0.; 276 254 cut = a*a - 4.*p*b; 277 if (cut > 0.) ecut = std::sqrt(cut);255 if (cut > 0.) { ecut = std::sqrt(cut); } 278 256 ecut = (ecut-a) / (p+p); 279 257 ecut2 = ecut; 280 //JMQ 290310 for avoiding unphysical increase below minimum (at ecut) 281 //ecut<0 means that there is no cut with energy axis, i.e. xs is set to 0 bellow minimum 282 // if (cut < 0.) ecut2 = ecut - 2.; 283 if (cut < 0.) ecut2 = ecut; 284 elab = K * FragmentA / ResidualA; 258 //JMQ 290310 for avoiding unphysical increase below minimum (at ecut) 259 // ecut<0 means that there is no cut with energy axis, i.e. xs is set 260 // to 0 bellow minimum 261 // if (cut < 0.) ecut2 = ecut - 2.; 262 if (cut < 0.) { ecut2 = ecut; } 263 elab = K * FragmentA /G4double(ResidualA); 285 264 sig = 0.; 286 265 if (elab <= ec) { //start for E<Ec 287 if (elab > ecut2) sig = (p*elab*elab+a*elab+b) * signor; 266 if (elab > ecut2) { sig = (p*elab*elab+a*elab+b) * signor; } 267 288 268 signor2 = (ec-elab-c) / w; 289 269 signor2 = 1. + std::exp(signor2); 290 sig = sig / signor2; 291 }//end for E<=Ec292 else 270 sig = sig / signor2; 271 } //end for E<=Ec 272 else{ //start for E>Ec 293 273 sig = (landa*elab+mu+nu/elab) * signor; 294 274 geom = 0.; 295 296 if (xnulam < flow || elab < etest) 297 {275 276 if (xnulam < flow || elab < etest) 277 { 298 278 if (sig <0.0) {sig=0.0;} 299 279 return sig; … … 303 283 geom = 31.416 * geom * geom; 304 284 sig = std::max(geom,sig); 305 285 306 286 } //end for E>Ec 307 return sig;} 308 309 310 311 // ************************** end of cross sections ******************************* 312 287 return sig; 288 } 289 -
trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4TritonEvaporationChannel.cc
r1340 r1347 25 25 // 26 26 // 27 // $Id: G4TritonEvaporationChannel.cc,v 1. 4 2006/06/29 20:10:45 gunterExp $28 // GEANT4 tag $Name: geant4-09-0 3-ref-09$27 // $Id: G4TritonEvaporationChannel.cc,v 1.5 2010/11/17 12:14:59 vnivanch Exp $ 28 // GEANT4 tag $Name: geant4-09-04-ref-00 $ 29 29 // 30 30 // Hadronic Process: Nuclear De-excitations 31 31 // by V. Lara (Nov. 1999) 32 32 // 33 // 17-11-2010 V.Ivanchenko moved constructor and destructor to source and cleanup 33 34 34 35 #include "G4TritonEvaporationChannel.hh" 35 36 37 G4TritonEvaporationChannel::G4TritonEvaporationChannel() 38 : G4EvaporationChannel(3,1,"triton",&theEvaporationProbability,&theCoulombBarrier) 39 {} 36 40 37 const G4TritonEvaporationChannel & G4TritonEvaporationChannel:: 38 operator=(const G4TritonEvaporationChannel & ) 39 { 40 throw G4HadronicException(__FILE__, __LINE__, "G4TritonEvaporationChannel::operator= meant to not be accessable"); 41 return *this; 42 } 41 G4TritonEvaporationChannel::~G4TritonEvaporationChannel() 42 {} 43 43 44 G4TritonEvaporationChannel::G4TritonEvaporationChannel(const G4TritonEvaporationChannel & ) : G4EvaporationChannel()45 {46 throw G4HadronicException(__FILE__, __LINE__, "G4TritonEvaporationChannel::CopyConstructor meant to not be accessable");47 }48 49 G4bool G4TritonEvaporationChannel::operator==(const G4TritonEvaporationChannel & right) const50 {51 return (this == (G4TritonEvaporationChannel *) &right);52 }53 54 G4bool G4TritonEvaporationChannel::operator!=(const G4TritonEvaporationChannel & right) const55 {56 return (this != (G4TritonEvaporationChannel *) &right);57 } -
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 -
trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4UnstableFragmentBreakUp.cc
r1340 r1347 25 25 // 26 26 // $Id: G4UnstableFragmentBreakUp.cc,v 1.8 2010/06/25 09:46:09 gunter Exp $ 27 // GEANT4 tag $Name: geant4-09-0 3-ref-09$27 // GEANT4 tag $Name: geant4-09-04-ref-00 $ 28 28 // 29 29 // ------------------------------------------------------------------- -
trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4VEvaporation.cc
r1340 r1347 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4VEvaporation.cc,v 1. 7 2010/04/27 14:00:40vnivanch Exp $27 // GEANT4 tag $Name: geant4-09-0 3-ref-09$26 // $Id: G4VEvaporation.cc,v 1.8 2010/10/29 17:35:03 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-04-ref-00 $ 28 28 // 29 29 // Hadronic Process: Nuclear De-excitations … … 33 33 #include "G4VEvaporation.hh" 34 34 35 G4VEvaporation::G4VEvaporation() 35 G4VEvaporation::G4VEvaporation():OPTxs(3),useSICB(false) 36 36 {} 37 37
Note: See TracChangeset
for help on using the changeset viewer.