91 using namespace G4InuclParticleNames;
92 using namespace G4InuclSpecialFunctions;
98 22.5, 22.0, 21.0, 21.0, 20.0,
101 20.6, 20.6, 18.6, 15.8, 13.5, 6.5,
104 6.65, 6.22, 6.27, 6.5, 6.7, 6.2, 6.25, 5.9, 6.1, 5.75,
107 6.46, 5.7, 6.28, 5.8, 6.15, 5.6, 5.8, 5.2, 5.8,
110 6.2 , 5.9 , 5.9, 6.0, 5.8, 5.7, 5.4, 5.4,
113 5.6, 6.1, 5.57, 6.3, 5.5, 5.8, 4.7, 6.2, 6.4, 6.2,
116 6.5, 6.2, 6.5, 5.3, 6.4, 5.7, 5.7, 6.2, 5.7,
119 6.3, 5.8, 6.7, 5.8, 6.6, 6.1, 4.3,
122 6.2, 3.8, 5.6, 4.0, 4.0, 4.2, 4.2, 3.5 };
126 0.6761, 0.677, 0.6788, 0.6803, 0.685,
128 0.6889, 0.6914, 0.6991, 0.7068, 0.725, 0.7391,
130 0.74, 0.741, 0.742, 0.743, 0.744, 0.7509, 0.752, 0.7531, 0.7543, 0.7548,
132 0.7557, 0.7566, 0.7576,
134 0.7587, 0.7597, 0.7608, 0.762, 0.7632, 0.7644, 0.7675, 0.7686, 0.7697, 0.7709,
136 0.7714, 0.7721, 0.7723, 0.7733, 0.7743, 0.7753, 0.7764,
138 0.7775, 0.7786, 0.7801, 0.781, 0.7821, 0.7831, 0.7842, 0.7852,
140 0.7864, 0.7875, 0.7880, 0.7887, 0.7889, 0.7899, 0.7909, 0.7919, 0.7930,
142 0.7941, 0.7953, 0.7965, 0.7977, 0.7987, 0.7989,
144 0.7997, 0.8075, 0.8097, 0.8119, 0.8143, 0.8164, 0.8174, 0.8274 };
152 theParaMaker(verboseLevel), QFinterp(XREP) {
153 parms.first.resize(6,0.);
154 parms.second.resize(6,0.);
171 G4cout <<
" >>> G4EquilibriumEvaporator::deExcite" <<
G4endl;
179 const G4double prob_cut_off = 1.0e-15;
180 const G4double Q1[6] = { 0.0, 0.0, 2.23, 8.49, 7.72, 28.3 };
181 const G4int AN[6] = { 1, 1, 2, 3, 3, 4 };
182 const G4int Q[6] = { 0, 1, 1, 1, 2, 2 };
183 const G4double G[6] = { 2.0, 2.0, 6.0, 6.0, 6.0, 4.0 };
185 const G4double fisssion_cut = 1000.0;
186 const G4double cut_off_energy = 0.1;
189 const G4int itry_max = 1000;
190 const G4int itry_global_max = 1000;
192 const G4int itry_gam_max = 100;
203 toTheNucleiSystemRestFrame.
setBullet(dummy);
208 if (explosion(
A,
Z,
EEXS)) {
210 theBigBanger.
deExcite(target, output);
217 if (
EEXS < cut_off_energy) {
226 G4double coul_coeff = (
A >= 100.0) ? 1.4 : 1.2;
231 G4bool fission_open =
true;
232 G4int itry_global = 0;
234 while (try_again && itry_global < itry_global_max) {
243 G4cout <<
" A " <<
A <<
" Z " <<
Z <<
" mass " << nuc_mass
247 if (explosion(
A,
Z,
EEXS)) {
249 G4cout <<
" big bang in eql step " << itry_global <<
G4endl;
257 if (
EEXS < cut_off_energy) {
259 G4cout <<
" no energy for evaporation in eql step " << itry_global
272 const std::vector<G4double>& AK = parms.first;
273 const std::vector<G4double>& CPA = parms.second;
278 for (i = 0; i < 6; i++) {
281 u[i] = parlev * A1[i];
285 if (goodRemnant(A1[i], Z1[i])) {
287 V[i] = coul_coeff *
Z * Q[i] * AK[i] / (1.0 +
EEXS / E0) /
289 TM[i] =
EEXS - QB - V[i] * A / A1[i];
291 G4cout <<
" i=" << i <<
" : QB " << QB <<
" u " << u[i]
292 <<
" V " << V[i] <<
" TM " << TM[i] <<
G4endl;
301 if (TM[0] > cut_off_energy) {
303 W[0] = BE *
G4cbrt(A1[0]*A1[0]) * G[0] * AL;
304 G4double TM1 = 2.0 * std::sqrt(u[0] * TM[0]) - ue;
306 if (TM1 > huge_num) TM1 = huge_num;
307 else if (TM1 < small) TM1 = small;
313 for (i = 1; i < 6; i++) {
315 if (TM[i] > cut_off_energy) {
316 W[i] = BE *
G4cbrt(A1[i]*A1[i]) * G[i] * (1.0 + CPA[i]);
317 G4double TM1 = 2.0 * std::sqrt(u[i] * TM[i]) - ue;
319 if (TM1 > huge_num) TM1 = huge_num;
320 else if (TM1 < small) TM1 = small;
329 if (A >= 100.0 && fission_open) {
332 G4double X = 0.019316 * X2 / (1.0 - 1.79 * X1 *
X1);
333 G4double EF = EEXS - getQF(X, X2, A,
Z, EEXS);
337 G4double TM1 = 2.0 * std::sqrt(AF * EF) - ue;
339 if (TM1 > huge_num) TM1 = huge_num;
340 else if (TM1 < small) TM1 = small;
342 W[6] = BF *
G4Exp(TM1);
343 if (W[6] > fisssion_cut*W[0]) W[6] = fisssion_cut*W[0];
351 G4cout <<
" Evaporation probabilities: sum = " << prob_sum
352 <<
"\n\t n " << W[0] <<
" p " << W[1] <<
" d " << W[2]
353 <<
" He3 " << W[3] <<
"\n\t t " << W[4] <<
" alpha " << W[5]
354 <<
" fission " << W[6] <<
G4endl;
359 if (prob_sum < prob_cut_off) {
361 G4double T00 = 1.0 / (std::sqrt(u1 / UCR0) - 1.25 / UCR0);
367 while (EEXS > cut_off_energy && try_again) {
374 FMAX = (T04*T04*T04*T04) *
G4Exp((EEXS - T04) / T00);
376 FMAX = EEXS*EEXS*EEXS*
EEXS;
380 G4cout <<
" T04 " << T04 <<
" FMAX (EEXS^4) " << FMAX <<
G4endl;
383 while (itry < itry_max) {
386 X1 = (S*S*S*S) *
G4Exp((EEXS - S) / T00);
391 if (itry == itry_max) {
409 <<
" pz " <<
PEX.
pz() <<
" E " <<
PEX.
e() << G4endl
410 <<
" evaporate px " << mom.
px() <<
" py " << mom.
py()
411 <<
" pz " << mom.
pz() <<
" E " << mom.
e() <<
G4endl;
426 if (itry_gam == itry_gam_max) try_again =
false;
435 for (i = 0; i < 7; i++) {
443 if (icase < 0)
continue;
447 G4cout <<
" particle/light-ion escape ("
448 << (icase==0 ?
"neutron" : icase==1 ?
"proton" :
449 icase==2 ?
"deuteron" : icase==3 ?
"He3" :
450 icase==4 ?
"triton" : icase==5 ?
"alpha" :
"ERROR!")
454 G4double uc = 2.0 * std::sqrt(u[icase] * TM[icase]);
461 while (itry1 < itry_max && bad) {
467 while (itry < itry_max && EPR < 0.0) {
470 S = 0.5 * (uc * uc - uu * uu) / u[icase];
471 EPR = TM[icase] - S * A / (A - 1.0) + V[icase];
475 G4cout <<
" EPR " << EPR <<
" V " << V[icase]
476 <<
" S " << S <<
" EEXS " << EEXS <<
G4endl;
479 if (EPR > 0.0 && S > V[icase] && S < EEXS) {
481 G4cout <<
" escape itry1 " << itry1 <<
" icase "
482 << icase <<
" S (MeV) " << S <<
G4endl;
487 G4int ptype = 2 - icase;
495 G4double pmod = std::sqrt((2.0 * mass + S) * S);
503 <<
" pz " <<
PEX.
pz() <<
" E " <<
PEX.
e() << G4endl
504 <<
" evaporate px " << mom.
px() <<
" py " << mom.
py()
505 <<
" pz " << mom.
pz() <<
" E " << mom.
e() <<
G4endl;
513 if (EEXS_new < 0.0)
continue;
532 G4cout <<
" nucleus A " << AN[icase] <<
" Z " << Q[icase]
533 <<
" escape icase " << icase <<
G4endl;
540 G4double pmod = std::sqrt((2.0 * mass + S) * S);
548 <<
" pz " <<
PEX.
pz() <<
" E " <<
PEX.
e() << G4endl
549 <<
" evaporate px " << mom.
px() <<
" py " << mom.
py()
550 <<
" pz " << mom.
pz() <<
" E " << mom.
e() <<
G4endl;
558 if (EEXS_new < 0.0)
continue;
568 AN[icase], Q[icase], 0.*
GeV,
580 if (itry1 == itry_max || bad) try_again =
false;
583 G4cout <<
" fission: A " << A <<
" Z " <<
Z <<
" eexs " << EEXS
584 <<
" Wn " << W[0] <<
" Wf " << W[6] <<
G4endl;
588 fission_output.
reset();
604 fission_open =
false;
612 if (itry_global == itry_global_max) {
614 G4cout <<
" ! itry_global " << itry_global_max <<
G4endl;
638 G4cout <<
" >>> G4EquilibriumEvaporator::explosion? ";
644 G4bool bigb = (!(a >= 12 && z >= 0 && z < 3*(a-
z)) &&
653 G4bool G4EquilibriumEvaporator::goodRemnant(
G4int a,
656 G4cout <<
" >>> G4EquilibriumEvaporator::goodRemnant(" << a <<
"," << z
657 <<
")? " << (a>1 && z>0 && a>
z) <<
G4endl;
660 return a > 1 && z > 0 && a >
z;
669 G4cout <<
" >>> G4EquilibriumEvaporator::getQF ";
678 if (x < XMIN || x > XMAX) {
680 G4double FX = (0.73 + (3.33 * X1 - 0.66) *
X1) * (X1*X1*X1);
682 QFF = G0 * FX *
G4cbrt(a*a);
687 if (QFF < 0.0) QFF = 0.0;
700 G4cout <<
" >>> G4EquilibriumEvaporator::getAF" <<
G4endl;
704 G4double AF = 1.285 * (1.0 - e / 1100.0);
706 if(AF < 1.06) AF = 1.06;
716 G4cout <<
" >>> G4EquilibriumEvaporator::getPARLEVDEN" <<
G4endl;
727 G4cout <<
" >>> G4EquilibriumEvaporator::getE0" <<
G4endl;
virtual void setVerboseLevel(G4int verbose=0)
virtual ~G4EquilibriumEvaporator()
void addOutgoingParticle(const G4InuclElementaryParticle &particle)
static G4double getParticleMass(G4int type)
void setBullet(const G4InuclParticle *bullet)
G4LorentzVector backToTheLab(const G4LorentzVector &mom) const
void getTargetData(const G4Fragment &target)
G4EquilibriumEvaporator()
G4GLOB_DLL std::ostream G4cout
G4double interpolate(const G4double x, const G4double(&yb)[nBins]) const
G4double getNucleiMass() const
void boostToLabFrame(const G4LorentzConvertor &convertor)
G4LorentzVector generateWithRandomAngles(G4double p, G4double mass=0.)
virtual void setVerboseLevel(G4int verbose)
G4double G4cbrt(G4double x)
G4double G4Log(G4double x)
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
const G4Fragment & makeFragment(G4LorentzVector mom, G4int A, G4int Z, G4double EX=0.)
const std::vector< G4InuclNuclei > & getOutgoingNuclei() const
const std::vector< G4InuclElementaryParticle > & getOutgoingParticles() const
void addOutgoingNucleus(const G4InuclNuclei &nuclei)
G4int numberOfFragments() const
virtual G4bool validateOutput(const G4Fragment &target, G4CollisionOutput &output)
void addRecoilFragment(const G4Fragment *aFragment)
const G4Fragment & getRecoilFragment(G4int index=0) const
virtual void deExcite(const G4Fragment &target, G4CollisionOutput &output)
void toTheTargetRestFrame()
void getParams(G4double Z, std::pair< std::vector< G4double >, std::vector< G4double > > &parms)
virtual void deExcite(const G4Fragment &target, G4CollisionOutput &output)
G4double bindingEnergy(G4int A, G4int Z)
virtual void deExcite(const G4Fragment &target, G4CollisionOutput &output)
void setTarget(const G4InuclParticle *target)