Geant4 11.1.1
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
G4EquilibriumEvaporator Class Reference

#include <G4EquilibriumEvaporator.hh>

+ Inheritance diagram for G4EquilibriumEvaporator:

Public Member Functions

 G4EquilibriumEvaporator ()
 
virtual ~G4EquilibriumEvaporator ()
 
virtual void setVerboseLevel (G4int verbose)
 
virtual void deExcite (const G4Fragment &target, G4CollisionOutput &output)
 
- Public Member Functions inherited from G4CascadeDeexciteBase
 G4CascadeDeexciteBase (const char *name)
 
virtual ~G4CascadeDeexciteBase ()
 
virtual void setVerboseLevel (G4int verbose=0)
 
- Public Member Functions inherited from G4VCascadeDeexcitation
 G4VCascadeDeexcitation (const G4String &name)
 
virtual ~G4VCascadeDeexcitation ()
 
virtual void collide (G4InuclParticle *bullet, G4InuclParticle *target, G4CollisionOutput &globalOutput)
 
virtual void deExcite (const G4Fragment &fragment, G4CollisionOutput &output)=0
 
- Public Member Functions inherited from G4VCascadeCollider
 G4VCascadeCollider (const G4String &name, G4int verbose=0)
 
virtual ~G4VCascadeCollider ()
 
virtual void collide (G4InuclParticle *bullet, G4InuclParticle *target, G4CollisionOutput &output)=0
 
virtual void setVerboseLevel (G4int verbose=0)
 

Additional Inherited Members

- Protected Member Functions inherited from G4CascadeDeexciteBase
virtual G4bool explosion (const G4Fragment &target) const
 
virtual G4bool explosion (G4int A, G4int Z, G4double excitation) const
 
virtual G4bool validateOutput (const G4Fragment &target, G4CollisionOutput &output)
 
virtual G4bool validateOutput (const G4Fragment &target, const std::vector< G4InuclElementaryParticle > &particles)
 
virtual G4bool validateOutput (const G4Fragment &target, const std::vector< G4InuclNuclei > &fragments)
 
void getTargetData (const G4Fragment &target)
 
const G4FragmentmakeFragment (G4LorentzVector mom, G4int A, G4int Z, G4double EX=0.)
 
const G4FragmentmakeFragment (G4int A, G4int Z, G4double EX=0.)
 
- Protected Member Functions inherited from G4VCascadeCollider
virtual void setName (const G4String &name)
 
- Protected Attributes inherited from G4CascadeDeexciteBase
G4CascadeCheckBalancebalance
 
G4int A
 
G4int Z
 
G4LorentzVector PEX
 
G4double EEXS
 
G4Fragment aFragment
 
- Protected Attributes inherited from G4VCascadeCollider
G4String theName
 
G4int verboseLevel
 

Detailed Description

Definition at line 57 of file G4EquilibriumEvaporator.hh.

Constructor & Destructor Documentation

◆ G4EquilibriumEvaporator()

G4EquilibriumEvaporator::G4EquilibriumEvaporator ( )

Definition at line 152 of file G4EquilibriumEvaporator.cc.

153 : G4CascadeDeexciteBase("G4EquilibriumEvaporator"),
154 theParaMaker(verboseLevel), QFinterp(XREP) {
155 parms.first.resize(6,0.);
156 parms.second.resize(6,0.);
157}

◆ ~G4EquilibriumEvaporator()

G4EquilibriumEvaporator::~G4EquilibriumEvaporator ( )
virtual

Definition at line 159 of file G4EquilibriumEvaporator.cc.

159{}

Member Function Documentation

◆ deExcite()

void G4EquilibriumEvaporator::deExcite ( const G4Fragment target,
G4CollisionOutput output 
)
virtual

Implements G4VCascadeDeexcitation.

Definition at line 170 of file G4EquilibriumEvaporator.cc.

171 {
172 if (verboseLevel) {
173 G4cout << " >>> G4EquilibriumEvaporator::deExcite" << G4endl;
174 }
175
176 if (verboseLevel>1) G4cout << " evaporating target: \n" << target << G4endl;
177
178 // Implementation of the equilibium evaporation according to Dostrovsky
179 // (Phys. Rev. 116 (1959) 683.
180 // Note that pairing energy which is of order 1 MeV for odd-even and even-odd
181 // nuclei is not included
182
183 // Element in array: 0 : neutron
184 // 1 : proton
185 // 2 : deuteron
186 // 3 : triton
187 // 4 : 3He
188 // 5 : alpha
189
190 const G4double huge_num = 50.0;
191 const G4double small = -50.0;
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 }; // binding energy
194 const G4int AN[6] = { 1, 1, 2, 3, 3, 4 }; // mass number
195 const G4int Q[6] = { 0, 1, 1, 1, 2, 2 }; // charge
196 const G4double G[6] = { 2.0, 2.0, 6.0, 6.0, 6.0, 4.0 };
197 const G4double BE = 0.0063;
198 const G4double fisssion_cut = 1000.0;
199 const G4double cut_off_energy = 0.1;
200
201 const G4double BF = 0.0242;
202 const G4int itry_max = 1000;
203 const G4int itry_global_max = 1000;
204 const G4double small_ekin = 1.0e-6;
205 const G4int itry_gam_max = 100;
206
207 G4double W[8];
208 G4double u[6]; // Level density for each particle type: (Atarg - Aemitted)*parlev
209 G4double V[6]; // Coulomb energy
210 G4double TM[6]; // Maximum possible kinetic energy of emitted particle
211 G4int A1[6]; // A of remnant nucleus
212 G4int Z1[6]; // Z of remnant nucleus
213
214 getTargetData(target);
215 if (verboseLevel > 3) G4cout << " after noeq: eexs " << EEXS << G4endl;
216
217 G4InuclElementaryParticle dummy(small_ekin, proton);
218 G4LorentzConvertor toTheNucleiSystemRestFrame;
219 //*** toTheNucleiSystemRestFrame.setVerbose(verboseLevel);
220 toTheNucleiSystemRestFrame.setBullet(dummy);
221
222 G4LorentzVector ppout;
223
224 // See if fragment should just be dispersed
225 if (explosion(A, Z, EEXS)) {
226 if (verboseLevel > 1) G4cout << " big bang in eql start " << G4endl;
227 theBigBanger.deExcite(target, output);
228
229 validateOutput(target, output); // Check energy conservation
230 return;
231 }
232
233 // If nucleus is in ground state, no evaporation
234 if (EEXS < cut_off_energy) {
235 if (verboseLevel > 1) G4cout << " no energy for evaporation" << G4endl;
236 output.addRecoilFragment(target);
237
238 validateOutput(target, output); // Check energy conservation
239 return;
240 }
241
242 // Initialize evaporation attempts
243 G4double coul_coeff = (A >= 100.0) ? 1.4 : 1.2;
244
245 G4LorentzVector pin = PEX; // Save original target for testing
246
247 G4bool try_again = true;
248 G4bool fission_open = true;
249 G4int itry_global = 0;
250
251 /* Loop checking 08.06.2015 MHK */
252 while (try_again && itry_global < itry_global_max) {
253 itry_global++;
254
255 // Set rest frame of current (recoiling) nucleus
256 toTheNucleiSystemRestFrame.setTarget(PEX);
257 toTheNucleiSystemRestFrame.toTheTargetRestFrame();
258
259 if (verboseLevel > 2) {
261 G4cout << " A " << A << " Z " << Z << " mass " << nuc_mass
262 << " EEXS " << EEXS << G4endl;
263 }
264
265 if (explosion(A, Z, EEXS)) { // big bang
266 if (verboseLevel > 2)
267 G4cout << " big bang in eql step " << itry_global << G4endl;
268
269 theBigBanger.deExcite(makeFragment(PEX,A,Z,EEXS), output);
270
271 validateOutput(target, output); // Check energy conservation
272 return;
273 }
274
275 if (EEXS < cut_off_energy) { // Evaporation not possible
276 if (verboseLevel > 2)
277 G4cout << " no energy for evaporation in eql step " << itry_global
278 << G4endl;
279
280 try_again = false;
281 break;
282 }
283
284 // Normal evaporation chain
285 G4double E0 = getE0(A);
286 G4double parlev = getPARLEVDEN(A, Z);
287 G4double u1 = parlev * A;
288
289 theParaMaker.getParams(Z, parms);
290 const std::vector<G4double>& AK = parms.first;
291 const std::vector<G4double>& CPA = parms.second;
292
293 G4double DM0 = bindingEnergy(A,Z);
294 G4int i(0);
295
296 for (i = 0; i < 6; i++) {
297 A1[i] = A - AN[i];
298 Z1[i] = Z - Q[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])) {
304 G4double QB = DM0 - bindingEnergy(A1[i],Z1[i]) - Q1[i];
305 V[i] = coul_coeff * Z * Q[i] * AK[i] / (1.0 + EEXS / E0) /
306 (G4cbrt(A1[i]) + G4cbrt(AN[i]));
307 TM[i] = EEXS - QB - V[i] * A / A1[i];
308 if (verboseLevel > 3) {
309 G4cout << " i=" << i << " : QB " << QB << " u " << u[i]
310 << " V " << V[i] << " TM " << TM[i] << G4endl;
311 }
312 }
313 }
314
315 G4double ue = 2.0 * std::sqrt(u1 * EEXS);
316 G4double prob_sum = 0.0;
317
318 W[0] = 0.0;
319 if (TM[0] > cut_off_energy) {
320 G4double AL = getAL(A);
321 G4double A13 = G4cbrt(A1[0]);
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
328 W[0] *= G4Exp(TM1);
329 prob_sum += W[0];
330 }
331
332 for (i = 1; i < 6; i++) {
333 W[i] = 0.0;
334 if (TM[i] > cut_off_energy) {
335 G4double A13 = G4cbrt(A1[i]);
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
342 W[i] *= G4Exp(TM1);
343 prob_sum += W[i];
344 }
345 }
346
347 // fisson part
348 W[6] = 0.0;
349 if (A >= 100.0 && fission_open) {
350 G4double X2 = Z * Z / A;
351 G4double X1 = 1.0 - 2.0 * Z / A;
352 G4double X = 0.019316 * X2 / (1.0 - 1.79 * X1 * X1);
353 G4double EF = EEXS - getQF(X, X2, A, Z, EEXS);
354
355 if (EF > 0.0) {
356 G4double AF = u1 * getAF(X, A, Z, EEXS);
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
362 W[6] = BF * G4Exp(TM1);
363 if (W[6] > fisssion_cut*W[0]) W[6] = fisssion_cut*W[0];
364
365 prob_sum += W[6];
366 }
367 }
368
369 // again time to decide what next
370 if (verboseLevel > 2) {
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
377 G4int icase = -1;
378
379 if (prob_sum < prob_cut_off) { // photon emission chain
380 G4double UCR0 = 2.5 + 150.0 / A;
381 G4double T00 = 1.0 / (std::sqrt(u1 / UCR0) - 1.25 / UCR0);
382
383 if (verboseLevel > 3)
384 G4cout << " UCR0 " << UCR0 << " T00 " << T00 << G4endl;
385
386 G4int itry_gam = 0;
387 while (EEXS > cut_off_energy && try_again) {
388 itry_gam++;
389 G4int itry = 0;
390 G4double T04 = 4.0 * T00;
391 G4double FMAX;
392
393 if (T04 < EEXS) {
394 FMAX = (T04*T04*T04*T04) * G4Exp((EEXS - T04) / T00);
395 } else {
396 FMAX = EEXS*EEXS*EEXS*EEXS;
397 };
398
399 if (verboseLevel > 3)
400 G4cout << " T04 " << T04 << " FMAX (EEXS^4) " << FMAX << G4endl;
401
402 G4double S(0), X1(0);
403 while (itry < itry_max) {
404 itry++;
405 S = EEXS * inuclRndm();
406 X1 = (S*S*S*S) * G4Exp((EEXS - S) / T00);
407
408 if (X1 > FMAX * inuclRndm()) break;
409 };
410
411 if (itry == itry_max) { // Maximum attempts exceeded
412 try_again = false;
413 break;
414 }
415
416 if (verboseLevel > 2) G4cout << " photon escape ?" << G4endl;
417
418 if (S < EEXS) { // Valid evaporate
419 S /= GeV; // Convert to Bertini units
420
421 G4double pmod = S;
423
424 // Push evaporated particle into current rest frame
425 mom = toTheNucleiSystemRestFrame.backToTheLab(mom);
426
427 if (verboseLevel > 3) {
428 G4cout << " nucleus px " << PEX.px() << " py " << PEX.py()
429 << " pz " << PEX.pz() << " E " << PEX.e() << G4endl
430 << " evaporate px " << mom.px() << " py " << mom.py()
431 << " pz " << mom.pz() << " E " << mom.e() << G4endl;
432 }
433
434 PEX -= mom; // Remaining four-momentum
435 EEXS -= S*GeV; // New excitation energy (in MeV)
436
437 // NOTE: In-situ construction will be optimized away (no copying)
440
441 if (verboseLevel > 3)
442 G4cout << output.getOutgoingParticles().back() << G4endl;
443
444 ppout += mom;
445 } else {
446 if (itry_gam == itry_gam_max) try_again = false;
447 }
448 } // while (EEXS > cut_off
449 try_again = false;
450 } else { // if (prob_sum < prob_cut_off)
451 G4double SL = prob_sum * inuclRndm();
452 if (verboseLevel > 3) G4cout << " random SL " << SL << G4endl;
453
454 G4double S1 = 0.0;
455 for (i = 0; i < 7; i++) { // Select evaporation scenario
456 S1 += W[i];
457 if (SL <= S1) {
458 icase = i;
459 break;
460 };
461 };
462
463 if (icase < 0) continue; // Failed to choose scenario, try again
464
465 if (icase < 6) { // particle or light nuclei escape
466 if (verboseLevel > 2) {
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!")
471 << ")?" << G4endl;
472 }
473
474 G4double Xmax = (std::sqrt(u[icase]*TM[icase] + 0.25) - 0.5)/u[icase];
475 G4int itry1 = 0;
476 G4bool bad = true;
477
478 while (itry1 < itry_max && bad) {
479 itry1++;
480 G4int itry = 0;
481 G4double S = 0.0;
482 G4double X = 0.0;
483 G4double Ptest = 0.0;
484
485 while (itry < itry_max) {
486 itry++;
487
488 // Sampling from eq. 17 of Dostrovsky
489 X = G4UniformRand()*TM[icase];
490 Ptest = (X/Xmax)*G4Exp(-2.*u[icase]*Xmax + 2.*std::sqrt(u[icase]*(TM[icase] - X)));
491 if (G4UniformRand() < Ptest) {
492 S = X + V[icase];
493 break;
494 }
495 }
496
497 if (S > V[icase] && S < EEXS) { // real escape
498
499 if (verboseLevel > 2)
500 G4cout << " escape itry1 " << itry1 << " icase "
501 << icase << " S (MeV) " << S << G4endl;
502
503 S /= GeV; // Convert to Bertini units
504
505 if (icase < 2) { // particle escape
506 G4int ptype = 2 - icase;
507 if (verboseLevel > 2)
508 G4cout << " particle " << ptype << " escape" << G4endl;
509
510 // generate particle momentum
511 G4double mass =
513
514 G4double pmod = std::sqrt((2.0 * mass + S) * S);
516
517 // Push evaporated particle into current rest frame
518 mom = toTheNucleiSystemRestFrame.backToTheLab(mom);
519
520 if (verboseLevel > 2) {
521 G4cout << " nucleus px " << PEX.px() << " py " << PEX.py()
522 << " pz " << PEX.pz() << " E " << PEX.e() << G4endl
523 << " evaporate px " << mom.px() << " py " << mom.py()
524 << " pz " << mom.pz() << " E " << mom.e() << G4endl;
525 }
526
527 // New excitation energy depends on residual nuclear state
528 G4double mass_new =
529 G4InuclNuclei::getNucleiMass(A1[icase],Z1[icase]);
530
531 G4double EEXS_new = ((PEX-mom).m() - mass_new)*GeV;
532 if (EEXS_new < 0.0) continue; // Sanity check for new nucleus
533
534 PEX -= mom; // Remaining four-momentum
535 EEXS = EEXS_new;
536
537 A = A1[icase];
538 Z = Z1[icase];
539
540 // NOTE: In-situ construction optimizes away (no copying)
543
544 if (verboseLevel > 3)
545 G4cout << output.getOutgoingParticles().back() << G4endl;
546
547 ppout += mom;
548 bad = false;
549 } else { // if (icase < 2)
550 if (verboseLevel > 2) {
551 G4cout << " nucleus A " << AN[icase] << " Z " << Q[icase]
552 << " escape icase " << icase << G4endl;
553 }
554
555 G4double mass =
556 G4InuclNuclei::getNucleiMass(AN[icase],Q[icase]);
557
558 // generate particle momentum
559 G4double pmod = std::sqrt((2.0 * mass + S) * S);
561
562 // Push evaporated particle into current rest frame
563 mom = toTheNucleiSystemRestFrame.backToTheLab(mom);
564
565 if (verboseLevel > 2) {
566 G4cout << " nucleus px " << PEX.px() << " py " << PEX.py()
567 << " pz " << PEX.pz() << " E " << PEX.e() << G4endl
568 << " evaporate px " << mom.px() << " py " << mom.py()
569 << " pz " << mom.pz() << " E " << mom.e() << G4endl;
570 }
571
572 // New excitation energy depends on residual nuclear state
573 G4double mass_new =
574 G4InuclNuclei::getNucleiMass(A1[icase],Z1[icase]);
575
576 G4double EEXS_new = ((PEX-mom).m() - mass_new)*GeV;
577 if (EEXS_new < 0.0) continue; // Sanity check for new nucleus
578
579 PEX -= mom; // Remaining four-momentum
580 EEXS = EEXS_new;
581
582 A = A1[icase];
583 Z = Z1[icase];
584
585 // NOTE: In-situ constructor optimizes away (no copying)
587 AN[icase], Q[icase], 0.*GeV,
589
590 if (verboseLevel > 3)
591 G4cout << output.getOutgoingNuclei().back() << G4endl;
592
593 ppout += mom;
594 bad = false;
595 } // if (icase < 2)
596 } // if (EPR > 0.0 ...
597 } // while (itry1 ...
598
599 if (itry1 == itry_max || bad) try_again = false;
600 } else { // if (icase < 6)
601 if (verboseLevel > 2) {
602 G4cout << " fission: A " << A << " Z " << Z << " eexs " << EEXS
603 << " Wn " << W[0] << " Wf " << W[6] << G4endl;
604 }
605
606 // Catch fission output separately for verification
607 fission_output.reset();
608 theFissioner.deExcite(makeFragment(A,Z,EEXS), fission_output);
609
610 if (fission_output.numberOfFragments() == 2) { // fission ok
611 if (verboseLevel > 2) G4cout << " fission done in eql" << G4endl;
612
613 // Move fission fragments to lab frame for processing
614 fission_output.boostToLabFrame(toTheNucleiSystemRestFrame);
615
616 // Now evaporate the fission fragments individually
617 this->deExcite(fission_output.getRecoilFragment(0), output);
618 this->deExcite(fission_output.getRecoilFragment(1), output);
619
620 validateOutput(target, output); // Check energy conservation
621 return;
622 } else { // fission forbidden now
623 fission_open = false;
624 }
625 } // End of fission case
626 } // if (prob_sum < prob_cut_off)
627 } // while (try_again
628
629 // this time it's final nuclei
630
631 if (itry_global == itry_global_max) {
632 if (verboseLevel > 3) {
633 G4cout << " ! itry_global " << itry_global_max << G4endl;
634 }
635 }
636
637 G4LorentzVector pnuc = pin - ppout;
638
639 // NOTE: In-situ constructor will be optimized away (no copying)
640 output.addOutgoingNucleus(G4InuclNuclei(pnuc, A, Z, 0.,
642
643 if (verboseLevel > 3) {
644 G4cout << " remaining nucleus \n" << output.getOutgoingNuclei().back()
645 << G4endl;
646 }
647
648
649 validateOutput(target, output); // Check energy conservation
650 return;
651}
G4double S(G4double temp)
#define A13
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:180
double G4double
Definition: G4Types.hh:83
bool G4bool
Definition: G4Types.hh:86
int G4int
Definition: G4Types.hh:85
#define G4endl
Definition: G4ios.hh:57
G4GLOB_DLL std::ostream G4cout
#define G4UniformRand()
Definition: Randomize.hh:52
virtual void deExcite(const G4Fragment &target, G4CollisionOutput &output)
Definition: G4BigBanger.cc:68
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)
Definition: G4Fissioner.cc:65
static G4double getParticleMass(G4int type)
G4double getNucleiMass() const
void getParams(G4double Z, std::pair< std::vector< G4double >, std::vector< G4double > > &parms)
Definition: paraMaker.cc:62
void setBullet(const G4InuclParticle *bullet)
G4LorentzVector backToTheLab(const G4LorentzVector &mom) const
void setTarget(const G4InuclParticle *target)
#define W
Definition: crc32.c:84
G4double bindingEnergy(G4int A, G4int Z)
G4LorentzVector generateWithRandomAngles(G4double p, G4double mass=0.)

Referenced by G4CascadeDeexcitation::deExcite(), deExcite(), and G4EvaporationInuclCollider::deExcite().

◆ setVerboseLevel()

void G4EquilibriumEvaporator::setVerboseLevel ( G4int  verbose)
virtual

Reimplemented from G4CascadeDeexciteBase.

Definition at line 161 of file G4EquilibriumEvaporator.cc.

161 {
163 theFissioner.setVerboseLevel(verbose);
164 theBigBanger.setVerboseLevel(verbose);
165}
virtual void setVerboseLevel(G4int verbose=0)

Referenced by G4CascadeDeexcitation::setVerboseLevel().


The documentation for this class was generated from the following files: