Geant4  10.03.p01
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
G4EquilibriumEvaporator Class Reference

#include <G4EquilibriumEvaporator.hh>

Inheritance diagram for G4EquilibriumEvaporator:
Collaboration 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 ()
 
- Public Member Functions inherited from G4VCascadeDeexcitation
 G4VCascadeDeexcitation (const G4String &name)
 
virtual ~G4VCascadeDeexcitation ()
 
virtual void collide (G4InuclParticle *bullet, G4InuclParticle *target, G4CollisionOutput &globalOutput)
 
- Public Member Functions inherited from G4VCascadeCollider
 G4VCascadeCollider (const G4String &name, G4int verbose=0)
 
virtual ~G4VCascadeCollider ()
 

Additional Inherited Members

- Protected Member Functions inherited from G4CascadeDeexciteBase
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 58 of file G4EquilibriumEvaporator.hh.

Constructor & Destructor Documentation

G4EquilibriumEvaporator::G4EquilibriumEvaporator ( )

Definition at line 153 of file G4EquilibriumEvaporator.cc.

154  : G4CascadeDeexciteBase("G4EquilibriumEvaporator"),
155  theParaMaker(verboseLevel), QFinterp(XREP) {
156  parms.first.resize(6,0.);
157  parms.second.resize(6,0.);
158 }
G4CascadeDeexciteBase(const char *name)
G4EquilibriumEvaporator::~G4EquilibriumEvaporator ( )
virtual

Definition at line 160 of file G4EquilibriumEvaporator.cc.

160 {}

Member Function Documentation

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

Implements G4VCascadeDeexcitation.

Definition at line 171 of file G4EquilibriumEvaporator.cc.

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

Here is the call graph for this function:

Here is the caller graph for this function:

void G4EquilibriumEvaporator::setVerboseLevel ( G4int  verbose)
virtual

Reimplemented from G4CascadeDeexciteBase.

Definition at line 162 of file G4EquilibriumEvaporator.cc.

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

Here is the call graph for this function:

Here is the caller graph for this function:


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