Changeset 1315 for trunk/source/processes/hadronic/models/pre_equilibrium
- Timestamp:
- Jun 18, 2010, 11:42:07 AM (14 years ago)
- Location:
- trunk/source/processes/hadronic/models/pre_equilibrium
- Files:
-
- 7 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/source/processes/hadronic/models/pre_equilibrium/History
r1228 r1315 14 14 * Please list in reverse chronological order (last date on top) 15 15 --------------------------------------------------------------- 16 17 09-April 2010 V.Ivanchenko hadr-pre-V09-03-04 18 --------------------------------------------------- 19 G4PreCompoundProton, G4PreCompoundDeuteron, G4PreCompoundTriton, 20 G4PreCompoundHe3, G4PreCompoundAlpha - (JMQ) return back to published 21 varian of OPT3 (Kalbach) parameterization of inverse 22 cross section 23 24 01-March 2010 V.Ivanchenko hadr-pre-V09-03-03 25 --------------------------------------------------- 26 G4VPreCompoundIon, G4VPreCompoundNucleon - removed dummy classes 27 28 25-February 2010 V.Ivanchenko hadr-pre-V09-03-02 29 --------------------------------------------------- 30 G4PreCompoundEmission - more protections for numerical computations are added 31 32 17-February 2010 V.Ivanchenko hadr-pre-V09-03-01 33 --------------------------------------------------- 34 G4PreCompoundEmission - (JMQ) return back value of Ef as it was in 9.3 35 36 19-January 2010 V.Ivanchenko hadr-pre-V09-03-00 37 --------------------------------------------------- 38 G4PreCompoundEmission - (JMQ) added protection against unphysical value of 39 parameter "an" when exitation energy is huge 40 to avoid numerical problem 41 - simplified algorithm of sampling 16 42 17 43 21-November 2009 V.Ivanchenko hadr-pre-V09-02-07 … … 29 55 12-November 2009 G.Cosmo hadr-pre-V09-02-05 30 56 --------------------------------------------------- 31 G4Pre compoundEmission - fixed compilation error on Windows.57 G4PreCompoundEmission - fixed compilation error on Windows. 32 58 33 59 12-November 2009 G.Folger hadr-pre-V09-02-04 -
trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/src/G4PreCompoundAlpha.cc
r1228 r1315 25 25 // 26 26 // 27 // $Id: G4PreCompoundAlpha.cc,v 1. 5 2009/02/13 18:57:32vnivanch Exp $28 // GEANT4 tag $Name: geant4-09-0 3$27 // $Id: G4PreCompoundAlpha.cc,v 1.6 2010/04/09 14:06:17 vnivanch Exp $ 28 // GEANT4 tag $Name: geant4-09-04-beta-cand-01 $ 29 29 // 30 30 // ------------------------------------------------------------------- … … 39 39 // Modified: 40 40 // 21.08.2008 J. M. Quesada add choice of options 41 // 10.02.2009 J. M. Quesada set default opt142 41 // 43 42 … … 243 242 ecut = (ecut-a) / (p+p); 244 243 ecut2 = ecut; 245 if (cut < 0.) ecut2 = ecut - 2.; 244 //JMQ 290310 for avoiding unphysical increase below minimum (at ecut) 245 //ecut<0 means that there is no cut with energy axis, i.e. xs is set to 0 bellow minimum 246 // if (cut < 0.) ecut2 = ecut - 2.; 247 if (cut < 0.) ecut2 = ecut; 246 248 elab = K * FragmentA / ResidualA; 247 249 sig = 0.; -
trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/src/G4PreCompoundDeuteron.cc
r1228 r1315 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4PreCompoundDeuteron.cc,v 1. 5 2009/02/13 18:57:32vnivanch Exp $27 // GEANT4 tag $Name: geant4-09-0 3$26 // $Id: G4PreCompoundDeuteron.cc,v 1.6 2010/04/09 14:06:17 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-04-beta-cand-01 $ 28 28 // 29 29 // ------------------------------------------------------------------- … … 38 38 // Modified: 39 39 // 21.08.2008 J. M. Quesada add choice of options 40 // 10.02.2009 J. M. Quesada set default opt141 40 // 42 41 … … 234 233 ecut = (ecut-a) / (p+p); 235 234 ecut2 = ecut; 236 if (cut < 0.) ecut2 = ecut - 2.; 235 //JMQ 290310 for avoiding unphysical increase below minimum (at ecut) 236 //ecut<0 means that there is no cut with energy axis, i.e. xs is set to 0 bellow minimum 237 // if (cut < 0.) ecut2 = ecut - 2.; 238 if (cut < 0.) ecut2 = ecut; 237 239 elab = K * FragmentA / ResidualA; 238 240 sig = 0.; -
trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/src/G4PreCompoundEmission.cc
r1228 r1315 25 25 // 26 26 // 27 // $Id: G4PreCompoundEmission.cc,v 1.2 2 2009/11/13 17:40:14 gcosmoExp $28 // GEANT4 tag $Name: geant4-09-0 3$27 // $Id: G4PreCompoundEmission.cc,v 1.28 2010/02/25 10:27:36 vnivanch Exp $ 28 // GEANT4 tag $Name: geant4-09-04-beta-cand-01 $ 29 29 // 30 30 // ------------------------------------------------------------------- … … 38 38 // 39 39 // Modified: 40 // 15.01.2010 J.M.Quesada added protection against unphysical values of parameter an 41 // 19.01.2010 V.Ivanchenko simplified computation of parameter an, sample cosTheta 42 // instead of theta; protect all calls to sqrt 40 43 // 41 44 … … 46 49 #include "G4HETCEmissionFactory.hh" 47 50 48 const G4PreCompoundEmission & G4PreCompoundEmission::operator=(const G4PreCompoundEmission &) 49 { 50 throw G4HadronicException(__FILE__, __LINE__, "G4PreCompoundEmission::operator= meant to not be accessable"); 51 const G4PreCompoundEmission & 52 G4PreCompoundEmission::operator=(const G4PreCompoundEmission &) 53 { 54 throw G4HadronicException(__FILE__, __LINE__, 55 "G4PreCompoundEmission::operator= meant to not be accessable"); 51 56 return *this; 52 57 } … … 166 171 // check that Excitation energy is >= 0 167 172 G4double anU = RestMomentum.m()-theFragment->GetRestNuclearMass(); 168 if (anU < 0.0) throw G4HadronicException(__FILE__, __LINE__, "G4PreCompoundModel::DeExcite: Excitation energy less than 0!"); 169 170 171 173 if (anU < 0.0) { 174 throw G4HadronicException(__FILE__, __LINE__, 175 "G4PreCompoundModel::DeExcite: Excitation energy less than 0!"); 176 } 177 172 178 // Update nucleus parameters: 173 179 // -------------------------- … … 204 210 } 205 211 206 207 G4 ThreeVector G4PreCompoundEmission::AngularDistribution(G4VPreCompoundFragment * theFragment,208 209 const G4double KineticEnergyOfEmittedFragment) const212 G4ThreeVector 213 G4PreCompoundEmission::AngularDistribution(G4VPreCompoundFragment * theFragment, 214 const G4Fragment& aFragment, 215 const G4double kinEnergyOfEmittedFrag) const 210 216 { 211 217 G4double p = aFragment.GetNumberOfParticles(); 212 218 G4double h = aFragment.GetNumberOfHoles(); 213 219 G4double U = aFragment.GetExcitationEnergy(); 214 215 // Kinetic Energy of emitted fragment 216 // G4double KineticEnergyOfEmittedFragment = theFragment->GetKineticEnergy(aFragment); 220 221 G4double ekin = std::max(0.0, kinEnergyOfEmittedFrag); 217 222 218 223 // Emission particle separation energy 219 224 G4double Bemission = theFragment->GetBindingEnergy(); 220 225 221 226 // Fermi energy 222 227 G4double Ef = G4PreCompoundParameters::GetAddress()->GetFermiEnergy(); … … 225 230 // G4EvaporationLevelDensityParameter theLDP; 226 231 // G4double g = (6.0/pi2)*aFragment.GetA()* 227 // theLDP.LevelDensityParameter(static_cast<G4int>(aFragment.GetA()),static_cast<G4int>(aFragment.GetZ()),U); 228 G4double g = (6.0/pi2)*aFragment.GetA() *229 G4PreCompoundParameters::GetAddress()->GetLevelDensity();232 233 G4double g = (6.0/pi2)*aFragment.GetA() 234 *G4PreCompoundParameters::GetAddress()->GetLevelDensity(); 230 235 231 236 // Average exciton energy relative to bottom of nuclear well … … 236 241 // G4double Uf = U - KineticEnergyOfEmittedFragment - Bemission; 237 242 238 239 240 243 G4double w_num = rho(p+1,h,g,Uf,Ef); 241 244 G4double w_den = rho(p,h,g,Uf,Ef); … … 251 254 252 255 253 G4double zeta = std::max(1.0,9.3/std::sqrt(KineticEnergyOfEmittedFragment/MeV)); 254 255 G4double an = 3.0*std::sqrt((ProjEnergy+Ef)*(KineticEnergyOfEmittedFragment+Bemission+Ef)); 256 if (aFragment.GetNumberOfExcitons() == 1) 257 { 258 an /= (zeta*2.0*aFragment.GetNumberOfExcitons()*Eav); 259 } 260 else 261 { 262 an /= (zeta*(aFragment.GetNumberOfExcitons()-1.0)*Eav); 263 } 256 // VI + JMQ 19/01/2010 update computation of the parameter an 257 // 258 G4double an = 0.0; 259 G4double Eeff = ekin + Bemission + Ef; 260 if(ekin > DBL_MIN && Eeff > DBL_MIN) { 261 262 G4double zeta = std::max(1.0,9.3/std::sqrt(ekin/MeV)); 263 264 an = 3.0*std::sqrt((ProjEnergy+Ef)*Eeff)/(zeta*Eav); 265 266 G4int ne = aFragment.GetNumberOfExcitons() - 1; 267 if ( ne > 1 ) { an /= (G4double)ne; } 264 268 265 266 // G4double expan=std::exp(an); 267 // thetaold = std::acos(std::log(expan-random*(expan-1.0/expan))/an); 268 // exp(an) becomes large for large an 269 // 1/exp(an) becomes large for large negative an 270 // take the large value out of the log depending on the sign of an to avoid FP exceptions 269 // protection of exponent 270 if ( an > 10. ) { an = 10.; } 271 } 272 273 // sample cosine of theta and not theta as in old versions 274 G4double random = G4UniformRand(); 275 G4double cost; 271 276 272 G4double random=G4UniformRand(); 273 G4double exp2an=0; 274 G4double theta; 275 if (an > 0.) { 276 if (an < 25.) { exp2an = std::exp(-2*an); } // we subtract from 1, exp(-50)~1e-21, so will not 277 // change numerical result 278 theta = std::acos(1+ std::log(1-random*(1-exp2an))/an); 279 } else if ( an < 0.) { 280 if ( an > -25.) { exp2an = std::exp(2*an); } // similar to above, except we compare to rndm*1 281 theta = std::acos(std::log(exp2an-random*(exp2an-1))/an - 1.); 282 } else { // an==0 now. 283 theta=std::acos(1.-2*random); 284 } 285 286 if (std::abs(an) < 50 ) 287 { 288 289 } 277 if(an < 0.1) { cost = 1. - 2*random; } 278 else { 279 G4double exp2an = std::exp(-2*an); 280 cost = 1. + std::log(1-random*(1-exp2an))/an; 281 if(cost > 1.) { cost = 1.; } 282 else if(cost < -1.) {cost = -1.; } 283 } 290 284 291 285 G4double phi = twopi*G4UniformRand(); 292 286 293 287 // Calculate the momentum magnitude of emitted fragment 294 G4double EmittedMass = theFragment->GetNuclearMass(); 295 G4double pmag = std::sqrt(KineticEnergyOfEmittedFragment*(KineticEnergyOfEmittedFragment+2.0*EmittedMass)); 296 297 298 G4double sinTheta = std::sin(theta); 299 // G4double cosTheta = std::sqrt(1.0-sinTheta*sinTheta); 300 G4double cosTheta = std::cos(theta); 301 302 G4ThreeVector momentum(pmag*std::cos(phi)*sinTheta,pmag*std::sin(phi)*sinTheta,pmag*cosTheta); 288 G4double pmag = std::sqrt(ekin*(ekin + 2.0*theFragment->GetNuclearMass())); 289 290 G4double sint = std::sqrt((1.0-cost)*(1.0+cost)); 291 292 G4ThreeVector momentum(pmag*std::cos(phi)*sint,pmag*std::sin(phi)*sint,pmag*cost); 303 293 // theta is the angle wrt the incident direction 304 294 momentum.rotateUz(theIncidentDirection); … … 309 299 G4double G4PreCompoundEmission::rho(const G4double p, const G4double h, const G4double g, 310 300 const G4double E, const G4double Ef) const 311 { 312 313 G4double Aph = (p*p + h*h + p - 3.0*h)/(4.0*g);314 G4double alpha = (p*p+h*h)/(2.0*g);315 316 if ( (E-alpha) < 0 ) return 0;301 { 302 // 25.02.2010 V.Ivanchenko added more protections 303 G4double Aph = (p*p + h*h + p - 3.0*h)/(4.0*g); 304 // G4double alpha = (p*p + h*h)/(2.0*g); 305 306 if ( E - Aph < 0.0) { return 0.0; } 317 307 318 308 G4double logConst = (p+h)*std::log(g) - logfactorial(p+h-1) - logfactorial(p) - logfactorial(h); 319 309 320 // initialise values using j=0310 // initialise values using j=0 321 311 322 312 G4double t1=1; 323 313 G4double t2=1; 324 G4double logt3=(p+h-1) * std::log(E-Aph); 325 G4double tot = std::exp( logt3 + logConst ); 326 327 // and now sum rest of terms 328 G4int j(1); 329 while ( (j <= h) && ((E - alpha - j*Ef) > 0.0) ) 330 { 331 t1 *= -1.; 332 t2 *= (h+1-j)/j; 333 logt3 = (p+h-1) * std::log( E - j*Ef - Aph) + logConst; 334 G4double t3 = std::exp(logt3); 335 tot += t1*t2*t3; 336 j++; 314 G4double logt3=(p+h-1) * std::log(E-Aph) + logConst; 315 const G4double logmax = 200.; 316 if(logt3 > logmax) { logt3 = logmax; } 317 G4double tot = std::exp( logt3 ); 318 319 // and now sum rest of terms 320 // 25.02.2010 V.Ivanchenko change while to for loop and cleanup 321 G4double Eeff = E - Aph; 322 for(G4int j=1; j<=h; ++j) 323 { 324 Eeff -= Ef; 325 if(Eeff < 0.0) { break; } 326 t1 *= -1.; 327 t2 *= (G4double)(h+1-j)/(G4double)j; 328 logt3 = (p+h-1) * std::log( Eeff) + logConst; 329 if(logt3 > logmax) { logt3 = logmax; } 330 tot += t1*t2*std::exp(logt3); 337 331 } 338 332 339 333 return tot; 340 334 } 341 342 343 335 344 336 G4double G4PreCompoundEmission::factorial(G4double a) const -
trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/src/G4PreCompoundHe3.cc
r1228 r1315 25 25 // 26 26 // 27 // $Id: G4PreCompoundHe3.cc,v 1. 5 2009/02/13 18:57:32vnivanch Exp $28 // GEANT4 tag $Name: geant4-09-0 3$27 // $Id: G4PreCompoundHe3.cc,v 1.6 2010/04/09 14:06:17 vnivanch Exp $ 28 // GEANT4 tag $Name: geant4-09-04-beta-cand-01 $ 29 29 // 30 30 // ------------------------------------------------------------------- … … 39 39 // Modified: 40 40 // 21.08.2008 J. M. Quesada add choice of options 41 // 10.02.2009 J. M. Quesada set default opt142 41 // 43 42 … … 242 241 ecut = (ecut-a) / (p+p); 243 242 ecut2 = ecut; 244 if (cut < 0.) ecut2 = ecut - 2.; 243 //JMQ 290310 for avoiding unphysical increase below minimum (at ecut) 244 //ecut<0 means that there is no cut with energy axis, i.e. xs is set to 0 bellow minimum 245 // if (cut < 0.) ecut2 = ecut - 2.; 246 if (cut < 0.) ecut2 = ecut; 245 247 elab = K * FragmentA / ResidualA; 246 248 sig = 0.; -
trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/src/G4PreCompoundProton.cc
r1228 r1315 25 25 // 26 26 // 27 // $Id: G4PreCompoundProton.cc,v 1. 4 2009/02/11 18:06:00vnivanch Exp $28 // GEANT4 tag $Name: geant4-09-0 3$27 // $Id: G4PreCompoundProton.cc,v 1.5 2010/04/09 14:06:17 vnivanch Exp $ 28 // GEANT4 tag $Name: geant4-09-04-beta-cand-01 $ 29 29 // 30 30 // ------------------------------------------------------------------- … … 68 68 //OPT=0 Dostrovski's parameterization 69 69 //OPT=1 Chatterjee's paramaterization 70 //OPT=2 Wellisch's parametarization70 //OPT=2,4 Wellisch's parametarization 71 71 //OPT=3 Kalbach's parameterization 72 72 // … … 224 224 // ** p from becchetti and greenlees (but modified with sub-barrier 225 225 // ** correction function and xp2 changed from -449) 226 // JMQ (june 2008) : refinement of proton cross section for light systems 227 // 226 228 227 G4double landa, landa0, landa1, mu, mu0, mu1,nu, nu0, nu1, nu2; 229 228 G4double p, p0, p1, p2; … … 281 280 ecut = (ecut-a) / (p+p); 282 281 ecut2 = ecut; 283 if (cut < 0.) ecut2 = ecut - 2.; 282 //JMQ 290310 for avoiding unphysical increase below minimum (at ecut) 283 //ecut<0 means that there is no cut with energy axis, i.e. xs is set to 0 bellow minimum 284 // if (cut < 0.) ecut2 = ecut - 2.; 285 if (cut < 0.) ecut2 = ecut; 284 286 elab = K * FragmentA / ResidualA; 285 287 sig = 0.; … … 290 292 signor2 = 1. + std::exp(signor2); 291 293 sig = sig / signor2; 292 // signor2 is empirical global corr'n at low elab for protons in PRECO, not enough for light targets 293 // refinement for proton cross section 294 if (ResidualZ<=26) 295 sig = sig*std::exp(-(a2*ResidualZ + b2)*(elab-(afit*ResidualZ+bfit)*ec)*(elab-(afit*ResidualZ+bfit)*ec)); 296 297 } //end for E<Ec 294 } //end for E<=Ec 298 295 else{ //start for E>Ec 299 296 sig = (landa*elab+mu+nu/elab) * signor; 300 301 // refinement for proton cross section302 if ( ResidualZ<=26 && elab <=(afit*ResidualZ+bfit)*ec)303 sig = sig*std::exp(-(a2*ResidualZ + b2)*(elab-(afit*ResidualZ+bfit)*ec)*(elab-(afit*ResidualZ+bfit)*ec));304 //305 306 297 geom = 0.; 307 298 -
trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/src/G4PreCompoundTriton.cc
r1228 r1315 25 25 // 26 26 // 27 // $Id: G4PreCompoundTriton.cc,v 1. 5 2009/02/13 18:57:32vnivanch Exp $28 // GEANT4 tag $Name: geant4-09-0 3$27 // $Id: G4PreCompoundTriton.cc,v 1.6 2010/04/09 14:06:17 vnivanch Exp $ 28 // GEANT4 tag $Name: geant4-09-04-beta-cand-01 $ 29 29 // 30 30 // ------------------------------------------------------------------- … … 39 39 // Modified: 40 40 // 21.08.2008 J. M. Quesada add choice of options 41 // 10.02.2009 J. M. Quesada set default opt142 41 // 43 42 … … 235 234 ecut = (ecut-a) / (p+p); 236 235 ecut2 = ecut; 237 if (cut < 0.) ecut2 = ecut - 2.; 236 //JMQ 290310 for avoiding unphysical increase below minimum (at ecut) 237 //ecut<0 means that there is no cut with energy axis, i.e. xs is set to 0 bellow minimum 238 // if (cut < 0.) ecut2 = ecut - 2.; 239 if (cut < 0.) ecut2 = ecut; 238 240 elab = K * FragmentA / ResidualA; 239 241 sig = 0.;
Note: See TracChangeset
for help on using the changeset viewer.