171 {
173 G4cout <<
" >>> G4EquilibriumEvaporator::deExcite" <<
G4endl;
174 }
175
177
178
179
180
181
182
183
184
185
186
187
188
189
192 const G4double prob_cut_off = 1.0e-15;
193 const G4double Q1[6] = { 0.0, 0.0, 2.23, 8.49, 7.72, 28.3 };
194 const G4int AN[6] = { 1, 1, 2, 3, 3, 4 };
195 const G4int Q[6] = { 0, 1, 1, 1, 2, 2 };
196 const G4double G[6] = { 2.0, 2.0, 6.0, 6.0, 6.0, 4.0 };
198 const G4double fisssion_cut = 1000.0;
199 const G4double cut_off_energy = 0.1;
200
202 const G4int itry_max = 1000;
203 const G4int itry_global_max = 1000;
205 const G4int itry_gam_max = 100;
206
213
216
219
220 toTheNucleiSystemRestFrame.
setBullet(dummy);
221
223
224
225 if (explosion(
A,
Z,
EEXS)) {
227 theBigBanger.
deExcite(target, output);
228
230 return;
231 }
232
233
234 if (
EEXS < cut_off_energy) {
237
239 return;
240 }
241
242
243 G4double coul_coeff = (
A >= 100.0) ? 1.4 : 1.2;
244
246
248 G4bool fission_open =
true;
249 G4int itry_global = 0;
250
251
252 while (try_again && itry_global < itry_global_max) {
253 itry_global++;
254
255
258
261 G4cout <<
" A " <<
A <<
" Z " <<
Z <<
" mass " << nuc_mass
263 }
264
265 if (explosion(
A,
Z,
EEXS)) {
267 G4cout <<
" big bang in eql step " << itry_global <<
G4endl;
268
270
272 return;
273 }
274
275 if (
EEXS < cut_off_energy) {
277 G4cout <<
" no energy for evaporation in eql step " << itry_global
279
280 try_again = false;
281 break;
282 }
283
284
288
290 const std::vector<G4double>& AK = parms.first;
291 const std::vector<G4double>& CPA = parms.second;
292
295
296 for (i = 0; i < 6; i++) {
299 u[i] = parlev * A1[i];
300 V[i] = 0.0;
301 TM[i] = -0.1;
302
303 if (goodRemnant(A1[i], Z1[i])) {
305 V[i] = coul_coeff *
Z * Q[i] * AK[i] / (1.0 +
EEXS / E0) /
307 TM[i] =
EEXS - QB - V[i] *
A / A1[i];
309 G4cout <<
" i=" << i <<
" : QB " << QB <<
" u " << u[i]
310 <<
" V " << V[i] <<
" TM " << TM[i] <<
G4endl;
311 }
312 }
313 }
314
317
319 if (TM[0] > cut_off_energy) {
322 W[0] = BE *
A13*
A13 * G[0] * AL;
323 G4double TM1 = 2.0 * std::sqrt(u[0] * TM[0]) - ue;
324
325 if (TM1 > huge_num) TM1 = huge_num;
326 else if (TM1 < small) TM1 = small;
327
330 }
331
332 for (i = 1; i < 6; i++) {
334 if (TM[i] > cut_off_energy) {
336 W[i] = BE *
A13*
A13 * G[i] * (1.0 + CPA[i]);
337 G4double TM1 = 2.0 * std::sqrt(u[i] * TM[i]) - ue;
338
339 if (TM1 > huge_num) TM1 = huge_num;
340 else if (TM1 < small) TM1 = small;
341
344 }
345 }
346
347
349 if (
A >= 100.0 && fission_open) {
352 G4double X = 0.019316 * X2 / (1.0 - 1.79 * X1 * X1);
354
355 if (EF > 0.0) {
357 G4double TM1 = 2.0 * std::sqrt(AF * EF) - ue;
358
359 if (TM1 > huge_num) TM1 = huge_num;
360 else if (TM1 < small) TM1 = small;
361
363 if (
W[6] > fisssion_cut*
W[0])
W[6] = fisssion_cut*
W[0];
364
366 }
367 }
368
369
371 G4cout <<
" Evaporation probabilities: sum = " << prob_sum
372 <<
"\n\t n " <<
W[0] <<
" p " <<
W[1] <<
" d " <<
W[2]
373 <<
" He3 " <<
W[3] <<
"\n\t t " <<
W[4] <<
" alpha " <<
W[5]
374 <<
" fission " <<
W[6] <<
G4endl;
375 }
376
378
379 if (prob_sum < prob_cut_off) {
381 G4double T00 = 1.0 / (std::sqrt(u1 / UCR0) - 1.25 / UCR0);
382
385
387 while (
EEXS > cut_off_energy && try_again) {
388 itry_gam++;
392
394 FMAX = (T04*T04*T04*T04) *
G4Exp((
EEXS - T04) / T00);
395 } else {
397 };
398
400 G4cout <<
" T04 " << T04 <<
" FMAX (EEXS^4) " << FMAX <<
G4endl;
401
403 while (itry < itry_max) {
404 itry++;
407
409 };
410
411 if (itry == itry_max) {
412 try_again = false;
413 break;
414 }
415
417
420
423
424
426
430 <<
" evaporate px " << mom.
px() <<
" py " << mom.
py()
431 <<
" pz " << mom.
pz() <<
" E " << mom.
e() <<
G4endl;
432 }
433
436
437
440
443
444 ppout += mom;
445 } else {
446 if (itry_gam == itry_gam_max) try_again = false;
447 }
448 }
449 try_again = false;
450 } else {
453
455 for (i = 0; i < 7; i++) {
457 if (SL <= S1) {
458 icase = i;
459 break;
460 };
461 };
462
463 if (icase < 0) continue;
464
465 if (icase < 6) {
467 G4cout <<
" particle/light-ion escape ("
468 << (icase==0 ? "neutron" : icase==1 ? "proton" :
469 icase==2 ? "deuteron" : icase==3 ? "He3" :
470 icase==4 ? "triton" : icase==5 ? "alpha" : "ERROR!")
472 }
473
474 G4double Xmax = (std::sqrt(u[icase]*TM[icase] + 0.25) - 0.5)/u[icase];
477
478 while (itry1 < itry_max && bad) {
479 itry1++;
484
485 while (itry < itry_max) {
486 itry++;
487
488
490 Ptest = (X/Xmax)*
G4Exp(-2.*u[icase]*Xmax + 2.*std::sqrt(u[icase]*(TM[icase] - X)));
493 break;
494 }
495 }
496
497 if (
S > V[icase] &&
S <
EEXS) {
498
500 G4cout <<
" escape itry1 " << itry1 <<
" icase "
501 << icase <<
" S (MeV) " <<
S <<
G4endl;
502
504
505 if (icase < 2) {
506 G4int ptype = 2 - icase;
509
510
513
514 G4double pmod = std::sqrt((2.0 * mass +
S) *
S);
516
517
519
523 <<
" evaporate px " << mom.
px() <<
" py " << mom.
py()
524 <<
" pz " << mom.
pz() <<
" E " << mom.
e() <<
G4endl;
525 }
526
527
530
531 G4double EEXS_new = ((
PEX-mom).m() - mass_new)*GeV;
532 if (EEXS_new < 0.0) continue;
533
536
539
540
543
546
547 ppout += mom;
548 bad = false;
549 } else {
551 G4cout <<
" nucleus A " << AN[icase] <<
" Z " << Q[icase]
552 <<
" escape icase " << icase <<
G4endl;
553 }
554
557
558
559 G4double pmod = std::sqrt((2.0 * mass +
S) *
S);
561
562
564
568 <<
" evaporate px " << mom.
px() <<
" py " << mom.
py()
569 <<
" pz " << mom.
pz() <<
" E " << mom.
e() <<
G4endl;
570 }
571
572
575
576 G4double EEXS_new = ((
PEX-mom).m() - mass_new)*GeV;
577 if (EEXS_new < 0.0) continue;
578
581
584
585
587 AN[icase], Q[icase], 0.*GeV,
589
592
593 ppout += mom;
594 bad = false;
595 }
596 }
597 }
598
599 if (itry1 == itry_max || bad) try_again = false;
600 } else {
602 G4cout <<
" fission: A " <<
A <<
" Z " <<
Z <<
" eexs " <<
EEXS
603 <<
" Wn " <<
W[0] <<
" Wf " <<
W[6] <<
G4endl;
604 }
605
606
607 fission_output.
reset();
609
612
613
615
616
619
621 return;
622 } else {
623 fission_open = false;
624 }
625 }
626 }
627 }
628
629
630
631 if (itry_global == itry_global_max) {
633 G4cout <<
" ! itry_global " << itry_global_max <<
G4endl;
634 }
635 }
636
638
639
642
646 }
647
648
650 return;
651}
G4double S(G4double temp)
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
G4GLOB_DLL std::ostream G4cout
virtual void deExcite(const G4Fragment &target, G4CollisionOutput &output)
void getTargetData(const G4Fragment &target)
const G4Fragment & makeFragment(G4LorentzVector mom, G4int A, G4int Z, G4double EX=0.)
virtual G4bool validateOutput(const G4Fragment &target, G4CollisionOutput &output)
void addRecoilFragment(const G4Fragment *aFragment)
const std::vector< G4InuclNuclei > & getOutgoingNuclei() const
void boostToLabFrame(const G4LorentzConvertor &convertor)
void addOutgoingParticle(const G4InuclElementaryParticle &particle)
const std::vector< G4InuclElementaryParticle > & getOutgoingParticles() const
const G4Fragment & getRecoilFragment(G4int index=0) const
void addOutgoingNucleus(const G4InuclNuclei &nuclei)
G4int numberOfFragments() const
virtual void deExcite(const G4Fragment &target, G4CollisionOutput &output)
virtual void deExcite(const G4Fragment &target, G4CollisionOutput &output)
static G4double getParticleMass(G4int type)
G4double getNucleiMass() const
void getParams(G4double Z, std::pair< std::vector< G4double >, std::vector< G4double > > &parms)
void toTheTargetRestFrame()
void setBullet(const G4InuclParticle *bullet)
G4LorentzVector backToTheLab(const G4LorentzVector &mom) const
void setTarget(const G4InuclParticle *target)
G4double bindingEnergy(G4int A, G4int Z)
G4LorentzVector generateWithRandomAngles(G4double p, G4double mass=0.)
G4double G4cbrt(G4double x)