95 {
97 G4cout <<
" >>> G4EquilibriumEvaporator::collide" <<
G4endl;
98 }
99
100
102 if (!nuclei_target) {
103 G4cerr <<
" EquilibriumEvaporator -> target is not nuclei " <<
G4endl;
104 return;
105 }
106
108
111
112
115 const G4double prob_cut_off = 1.0e-15;
116 const G4double Q1[6] = { 0.0, 0.0, 2.23, 8.49, 7.72, 28.3 };
117 const G4int AN[6] = { 1, 1, 2, 3, 3, 4 };
118 const G4int Q[6] = { 0, 1, 1, 1, 2, 2 };
119 const G4double G[6] = { 2.0, 2.0, 6.0, 6.0, 6.0, 4.0 };
121 const G4double fisssion_cut = 1000.0;
122 const G4double cut_off_energy = 0.1;
123
125 const G4int itry_max = 1000;
126 const G4int itry_global_max = 1000;
128 const G4int itry_gam_max = 100;
129
132
137
139
142
143 toTheNucleiSystemRestFrame.
setBullet(dummy);
144
146
147
148 if (explosion(A, Z, EEXS)) {
150 theBigBanger.
collide(0, target, output);
151
153 return;
154 }
155
156
157 if (EEXS < cut_off_energy) {
160
162 return;
163 }
164
165
166 G4double coul_coeff = (A >= 100.0) ? 1.4 : 1.2;
167
169
171 G4bool fission_open =
true;
172 G4int itry_global = 0;
173
174 while (try_again && itry_global < itry_global_max) {
175 itry_global++;
176
177
178 toTheNucleiSystemRestFrame.
setTarget(PEX);
180
183 G4cout <<
" A " << A <<
" Z " << Z <<
" mass " << nuc_mass
184 <<
" EEXS " << EEXS <<
G4endl;
185 }
186
187 if (explosion(A, Z, EEXS)) {
189 G4cout <<
" big bang in eql step " << itry_global <<
G4endl;
190
193
195 return;
196 }
197
198 if (EEXS < cut_off_energy) {
200 G4cout <<
" no energy for evaporation in eql step " << itry_global
202
203 try_again = false;
204 break;
205 }
206
207
209 G4double parlev = getPARLEVDEN(A, Z);
211
213 const std::vector<G4double>& AK = parms.first;
214 const std::vector<G4double>& CPA = parms.second;
215
218
219 for (i = 0; i < 6; i++) {
220 A1[i] = A - AN[i];
221 Z1[i] = Z - Q[i];
222 u[i] = parlev * A1[i];
223 V[i] = 0.0;
224 TM[i] = -0.1;
225
226 if (goodRemnant(A1[i], Z1[i])) {
228 V[i] = coul_coeff * Z * Q[i] * AK[i] / (1.0 + EEXS / E0) /
230 TM[i] = EEXS - QB - V[i] * A / A1[i];
232 G4cout <<
" i=" << i <<
" : QB " << QB <<
" u " << u[i]
233 <<
" V " << V[i] <<
" TM " << TM[i] <<
G4endl;
234 }
235 }
236 }
237
238 G4double ue = 2.0 * std::sqrt(u1 * EEXS);
240
241 W[0] = 0.0;
242 if (TM[0] > cut_off_energy) {
244 W[0] = BE *
G4cbrt(A1[0]*A1[0]) * G[0] * AL;
245 G4double TM1 = 2.0 * std::sqrt(u[0] * TM[0]) - ue;
246
247 if (TM1 > huge_num) TM1 = huge_num;
248 else if (TM1 < small) TM1 = small;
249
250 W[0] *= std::exp(TM1);
251 prob_sum += W[0];
252 }
253
254 for (i = 1; i < 6; i++) {
255 W[i] = 0.0;
256 if (TM[i] > cut_off_energy) {
257 W[i] = BE *
G4cbrt(A1[i]*A1[i]) * G[i] * (1.0 + CPA[i]);
258 G4double TM1 = 2.0 * std::sqrt(u[i] * TM[i]) - ue;
259
260 if (TM1 > huge_num) TM1 = huge_num;
261 else if (TM1 < small) TM1 = small;
262
263 W[i] *= std::exp(TM1);
264 prob_sum += W[i];
265 }
266 }
267
268
269 W[6] = 0.0;
270 if (A >= 100.0 && fission_open) {
273 G4double X = 0.019316 * X2 / (1.0 - 1.79 * X1 * X1);
274 G4double EF = EEXS - getQF(X, X2, A, Z, EEXS);
275
276 if (EF > 0.0) {
277 G4double AF = u1 * getAF(X, A, Z, EEXS);
278 G4double TM1 = 2.0 * std::sqrt(AF * EF) - ue;
279
280 if (TM1 > huge_num) TM1 = huge_num;
281 else if (TM1 < small) TM1 = small;
282
283 W[6] = BF * std::exp(TM1);
284 if (W[6] > fisssion_cut*W[0]) W[6] = fisssion_cut*W[0];
285
286 prob_sum += W[6];
287 }
288 }
289
290
292 G4cout <<
" Evaporation probabilities: sum = " << prob_sum
293 << "\n\t n " << W[0] << " p " << W[1] << " d " << W[2]
294 << " He3 " << W[3] << "\n\t t " << W[4] << " alpha " << W[5]
295 <<
" fission " << W[6] <<
G4endl;
296 }
297
299
300 if (prob_sum < prob_cut_off) {
302 G4double T00 = 1.0 / (std::sqrt(u1 / UCR0) - 1.25 / UCR0);
303
306
308 while (EEXS > cut_off_energy && try_again) {
309 itry_gam++;
313
314 if (T04 < EEXS) {
315 FMAX = (T04*T04*T04*T04) * std::exp((EEXS - T04) / T00);
316 } else {
317 FMAX = EEXS*EEXS*EEXS*EEXS;
318 };
319
321 G4cout <<
" T04 " << T04 <<
" FMAX (EEXS^4) " << FMAX <<
G4endl;
322
324 while (itry < itry_max) {
325 itry++;
327 X1 = (S*S*S*S) * std::exp((EEXS - S) / T00);
328
330 };
331
332 if (itry == itry_max) {
333 try_again = false;
334 break;
335 }
336
338
339 if (S < EEXS) {
340 S /= GeV;
341
344
345
347
349 G4cout <<
" nucleus px " << PEX.
px() <<
" py " << PEX.
py()
350 <<
" pz " << PEX.
pz() <<
" E " << PEX.
e() <<
G4endl
351 <<
" evaporate px " << mom.
px() <<
" py " << mom.
py()
352 <<
" pz " << mom.
pz() <<
" E " << mom.
e() <<
G4endl;
353 }
354
355 PEX -= mom;
356 EEXS -= S*GeV;
357
358
361
364
365 ppout += mom;
366 } else {
367 if (itry_gam == itry_gam_max) try_again = false;
368 }
369 }
370 try_again = false;
371 } else {
374
376 for (i = 0; i < 7; i++) {
377 S1 += W[i];
378 if (SL <= S1) {
379 icase = i;
380 break;
381 };
382 };
383
384 if (icase < 0) continue;
385
386 if (icase < 6) {
388 G4cout <<
" particle/light-ion escape ("
389 << (icase==0 ? "neutron" : icase==1 ? "proton" :
390 icase==2 ? "deuteron" : icase==3 ? "He3" :
391 icase==4 ? "triton" : icase==5 ? "alpha" : "ERROR!")
393 }
394
395 G4double uc = 2.0 * std::sqrt(u[icase] * TM[icase]);
396 G4double ur = (uc > huge_num ? std::exp(huge_num) : std::exp(uc));
401
402 while (itry1 < itry_max && bad) {
403 itry1++;
407
408 while (itry < itry_max && EPR < 0.0) {
409 itry++;
411 S = 0.5 * (uc * uc - uu * uu) / u[icase];
412 EPR = TM[icase] - S * A / (A - 1.0) + V[icase];
413 };
414
416 G4cout <<
" EPR " << EPR <<
" V " << V[icase]
417 <<
" S " << S <<
" EEXS " << EEXS <<
G4endl;
418 }
419
420 if (EPR > 0.0 && S > V[icase] && S < EEXS) {
422 G4cout <<
" escape itry1 " << itry1 <<
" icase "
423 << icase <<
" S (MeV) " << S <<
G4endl;
424
425 S /= GeV;
426
427 if (icase < 2) {
428 G4int ptype = 2 - icase;
431
432
435
436 G4double pmod = std::sqrt((2.0 * mass + S) * S);
438
439
441
443 G4cout <<
" nucleus px " << PEX.
px() <<
" py " << PEX.
py()
444 <<
" pz " << PEX.
pz() <<
" E " << PEX.
e() <<
G4endl
445 <<
" evaporate px " << mom.
px() <<
" py " << mom.
py()
446 <<
" pz " << mom.
pz() <<
" E " << mom.
e() <<
G4endl;
447 }
448
449
452
453 G4double EEXS_new = ((PEX-mom).m() - mass_new)*GeV;
454 if (EEXS_new < 0.0) continue;
455
456 PEX -= mom;
457 EEXS = EEXS_new;
458
459 A = A1[icase];
460 Z = Z1[icase];
461
462
465
468
469 ppout += mom;
470 bad = false;
471 } else {
473 G4cout <<
" nucleus A " << AN[icase] <<
" Z " << Q[icase]
474 <<
" escape icase " << icase <<
G4endl;
475 }
476
479
480
481 G4double pmod = std::sqrt((2.0 * mass + S) * S);
483
484
486
488 G4cout <<
" nucleus px " << PEX.
px() <<
" py " << PEX.
py()
489 <<
" pz " << PEX.
pz() <<
" E " << PEX.
e() <<
G4endl
490 <<
" evaporate px " << mom.
px() <<
" py " << mom.
py()
491 <<
" pz " << mom.
pz() <<
" E " << mom.
e() <<
G4endl;
492 }
493
494
497
498 G4double EEXS_new = ((PEX-mom).m() - mass_new)*GeV;
499 if (EEXS_new < 0.0) continue;
500
501 PEX -= mom;
502 EEXS = EEXS_new;
503
504 A = A1[icase];
505 Z = Z1[icase];
506
507
509 AN[icase], Q[icase], 0.*GeV,
511
514
515 ppout += mom;
516 bad = false;
517 }
518 }
519 }
520
521 if (itry1 == itry_max || bad) try_again = false;
522 } else {
524
526 G4cout <<
" fission: A " << A <<
" Z " << Z <<
" eexs " << EEXS
527 <<
" Wn " << W[0] <<
" Wf " << W[6] <<
G4endl;
528 }
529
530
531 fission_output.
reset();
533
535 if (nuclea.size() == 2) {
537
538
540
541
544
545 this->
collide(0, &nuclea[0], output);
546 this->
collide(0, &nuclea[1], output);
547
550 return;
551 } else {
552 fission_open = false;
553 }
554 }
555 }
556 }
557
558
559
560 if (itry_global == itry_global_max) {
562 G4cout <<
" ! itry_global " << itry_global_max <<
G4endl;
563 }
564 }
565
567
568
571
575 }
576
577
579 return;
580}
G4DLLIMPORT std::ostream G4cerr
G4DLLIMPORT std::ostream G4cout
void collide(G4InuclParticle *bullet, G4InuclParticle *target, G4CollisionOutput &output)
virtual void setVerboseLevel(G4int verbose=0)
virtual void setConservationChecks(G4bool doBalance=true)
G4bool doConservationChecks
virtual G4bool validateOutput(G4InuclParticle *bullet, G4InuclParticle *target, G4CollisionOutput &output)
const std::vector< G4InuclNuclei > & getOutgoingNuclei() const
void boostToLabFrame(const G4LorentzConvertor &convertor)
void addOutgoingParticle(const G4InuclElementaryParticle &particle)
const std::vector< G4InuclElementaryParticle > & getOutgoingParticles() const
void addOutgoingNucleus(const G4InuclNuclei &nuclei)
void collide(G4InuclParticle *bullet, G4InuclParticle *target, G4CollisionOutput &output)
void collide(G4InuclParticle *bullet, G4InuclParticle *target, G4CollisionOutput &output)
static G4double getParticleMass(G4int type)
G4double getNucleiMass() const
G4double getExitationEnergy() const
G4LorentzVector getMomentum() const
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.)
void paraMaker(G4double Z, std::pair< std::vector< G4double >, std::vector< G4double > > &parms)
G4double G4cbrt(G4double x)