Geant4 9.6.0
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 ()
 
void collide (G4InuclParticle *bullet, G4InuclParticle *target, G4CollisionOutput &output)
 
- Public Member Functions inherited from G4CascadeColliderBase
 G4CascadeColliderBase (const char *name, G4int verbose=0)
 
virtual ~G4CascadeColliderBase ()
 
virtual void rescatter (G4InuclParticle *, G4KineticTrackVector *, G4V3DNucleus *, G4CollisionOutput &)
 
virtual void setVerboseLevel (G4int verbose=0)
 
virtual void setConservationChecks (G4bool doBalance=true)
 
- Public Member Functions inherited from G4VCascadeCollider
 G4VCascadeCollider (const char *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 G4CascadeColliderBase
virtual G4bool useEPCollider (G4InuclParticle *bullet, G4InuclParticle *target) const
 
virtual G4bool explosion (G4InuclNuclei *target) const
 
virtual G4bool explosion (G4Fragment *target) const
 
virtual G4bool explosion (G4int A, G4int Z, G4double excitation) const
 
virtual G4bool inelasticInteractionPossible (G4InuclParticle *bullet, G4InuclParticle *target, G4double ekin) const
 
virtual G4bool validateOutput (G4InuclParticle *bullet, G4InuclParticle *target, G4CollisionOutput &output)
 
virtual G4bool validateOutput (G4InuclParticle *bullet, G4InuclParticle *target, const std::vector< G4InuclElementaryParticle > &particles)
 
virtual G4bool validateOutput (G4InuclParticle *bullet, G4InuclParticle *target, const std::vector< G4InuclNuclei > &fragments)
 
- Protected Member Functions inherited from G4VCascadeCollider
virtual void setName (const char *name)
 
- Protected Attributes inherited from G4CascadeColliderBase
G4InteractionCase interCase
 
G4bool doConservationChecks
 
G4CascadeCheckBalancebalance
 
- Protected Attributes inherited from G4VCascadeCollider
const char * theName
 
G4int verboseLevel
 

Detailed Description

Definition at line 49 of file G4EquilibriumEvaporator.hh.

Constructor & Destructor Documentation

◆ G4EquilibriumEvaporator()

G4EquilibriumEvaporator::G4EquilibriumEvaporator ( )

Definition at line 84 of file G4EquilibriumEvaporator.cc.

85 : G4CascadeColliderBase("G4EquilibriumEvaporator") {
86 parms.first.resize(6,0.);
87 parms.second.resize(6,0.);
88}

◆ ~G4EquilibriumEvaporator()

G4EquilibriumEvaporator::~G4EquilibriumEvaporator ( )
virtual

Definition at line 90 of file G4EquilibriumEvaporator.cc.

90{}

Member Function Documentation

◆ collide()

void G4EquilibriumEvaporator::collide ( G4InuclParticle bullet,
G4InuclParticle target,
G4CollisionOutput output 
)
virtual

Implements G4VCascadeCollider.

Definition at line 93 of file G4EquilibriumEvaporator.cc.

95 {
96 if (verboseLevel) {
97 G4cout << " >>> G4EquilibriumEvaporator::collide" << G4endl;
98 }
99
100 // Sanity check
101 G4InuclNuclei* nuclei_target = dynamic_cast<G4InuclNuclei*>(target);
102 if (!nuclei_target) {
103 G4cerr << " EquilibriumEvaporator -> target is not nuclei " << G4endl;
104 return;
105 }
106
107 if (verboseLevel>1) G4cout << " evaporating target: \n" << *target << G4endl;
108
109 theFissioner.setVerboseLevel(verboseLevel);
110 theBigBanger.setVerboseLevel(verboseLevel);
111
112 // simple implementation of the equilibium evaporation a la Dostrowski
113 const G4double huge_num = 50.0;
114 const G4double small = -50.0;
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 };
120 const G4double BE = 0.0063;
121 const G4double fisssion_cut = 1000.0;
122 const G4double cut_off_energy = 0.1;
123
124 const G4double BF = 0.0242;
125 const G4int itry_max = 1000;
126 const G4int itry_global_max = 1000;
127 const G4double small_ekin = 1.0e-6;
128 const G4int itry_gam_max = 100;
129
130 G4double W[8], u[6], V[6], TM[6];
131 G4int A1[6], Z1[6];
132
133 G4int A = nuclei_target->getA();
134 G4int Z = nuclei_target->getZ();
135 G4LorentzVector PEX = nuclei_target->getMomentum();
136 G4double EEXS = nuclei_target->getExitationEnergy();
137
138 if (verboseLevel > 3) G4cout << " after noeq: eexs " << EEXS << G4endl;
139
140 G4InuclElementaryParticle dummy(small_ekin, proton);
141 G4LorentzConvertor toTheNucleiSystemRestFrame;
142 //*** toTheNucleiSystemRestFrame.setVerbose(verboseLevel);
143 toTheNucleiSystemRestFrame.setBullet(dummy);
144
145 G4LorentzVector ppout;
146
147 // See if fragment should just be dispersed
148 if (explosion(A, Z, EEXS)) {
149 if (verboseLevel > 1) G4cout << " big bang in eql start " << G4endl;
150 theBigBanger.collide(0, target, output);
151
152 validateOutput(0, target, output); // Check energy conservation
153 return;
154 }
155
156 // If nucleus is in ground state, no evaporation
157 if (EEXS < cut_off_energy) {
158 if (verboseLevel > 1) G4cout << " no energy for evaporation" << G4endl;
159 output.addOutgoingNucleus(*nuclei_target);
160
161 validateOutput(0, target, output); // Check energy conservation
162 return;
163 }
164
165 // Initialize evaporation attempts
166 G4double coul_coeff = (A >= 100.0) ? 1.4 : 1.2;
167
168 G4LorentzVector pin = PEX; // Save original target for testing
169
170 G4bool try_again = true;
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 // Set rest frame of current (recoiling) nucleus
178 toTheNucleiSystemRestFrame.setTarget(PEX);
179 toTheNucleiSystemRestFrame.toTheTargetRestFrame();
180
181 if (verboseLevel > 2) {
182 G4double nuc_mass = G4InuclNuclei::getNucleiMass(A, Z, EEXS);
183 G4cout << " A " << A << " Z " << Z << " mass " << nuc_mass
184 << " EEXS " << EEXS << G4endl;
185 }
186
187 if (explosion(A, Z, EEXS)) { // big bang
188 if (verboseLevel > 2)
189 G4cout << " big bang in eql step " << itry_global << G4endl;
190
192 theBigBanger.collide(0, &nuclei, output);
193
194 validateOutput(0, target, output); // Check energy conservation
195 return;
196 }
197
198 if (EEXS < cut_off_energy) { // Evaporation not possible
199 if (verboseLevel > 2)
200 G4cout << " no energy for evaporation in eql step " << itry_global
201 << G4endl;
202
203 try_again = false;
204 break;
205 }
206
207 // Normal evaporation chain
208 G4double E0 = getE0(A);
209 G4double parlev = getPARLEVDEN(A, Z);
210 G4double u1 = parlev * A;
211
212 paraMaker(Z, parms);
213 const std::vector<G4double>& AK = parms.first;
214 const std::vector<G4double>& CPA = parms.second;
215
216 G4double DM0 = bindingEnergy(A,Z);
217 G4int i(0);
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])) {
227 G4double QB = DM0 - bindingEnergy(A1[i],Z1[i]) - Q1[i];
228 V[i] = coul_coeff * Z * Q[i] * AK[i] / (1.0 + EEXS / E0) /
229 (G4cbrt(A1[i]) + G4cbrt(AN[i]));
230 TM[i] = EEXS - QB - V[i] * A / A1[i];
231 if (verboseLevel > 3) {
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);
239 G4double prob_sum = 0.0;
240
241 W[0] = 0.0;
242 if (TM[0] > cut_off_energy) {
243 G4double AL = getAL(A);
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 // fisson part
269 W[6] = 0.0;
270 if (A >= 100.0 && fission_open) {
271 G4double X2 = Z * Z / A;
272 G4double X1 = 1.0 - 2.0 * Z / A;
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 // again time to decide what next
291 if (verboseLevel > 2) {
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
298 G4int icase = -1;
299
300 if (prob_sum < prob_cut_off) { // photon emission chain
301 G4double UCR0 = 2.5 + 150.0 / A;
302 G4double T00 = 1.0 / (std::sqrt(u1 / UCR0) - 1.25 / UCR0);
303
304 if (verboseLevel > 3)
305 G4cout << " UCR0 " << UCR0 << " T00 " << T00 << G4endl;
306
307 G4int itry_gam = 0;
308 while (EEXS > cut_off_energy && try_again) {
309 itry_gam++;
310 G4int itry = 0;
311 G4double T04 = 4.0 * T00;
312 G4double FMAX;
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
320 if (verboseLevel > 3)
321 G4cout << " T04 " << T04 << " FMAX (EEXS^4) " << FMAX << G4endl;
322
323 G4double S(0), X1(0);
324 while (itry < itry_max) {
325 itry++;
326 S = EEXS * inuclRndm();
327 X1 = (S*S*S*S) * std::exp((EEXS - S) / T00);
328
329 if (X1 > FMAX * inuclRndm()) break;
330 };
331
332 if (itry == itry_max) { // Maximum attempts exceeded
333 try_again = false;
334 break;
335 }
336
337 if (verboseLevel > 2) G4cout << " photon escape ?" << G4endl;
338
339 if (S < EEXS) { // Valid evaporate
340 S /= GeV; // Convert to Bertini units
341
342 G4double pmod = S;
344
345 // Push evaporated particle into current rest frame
346 mom = toTheNucleiSystemRestFrame.backToTheLab(mom);
347
348 if (verboseLevel > 3) {
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; // Remaining four-momentum
356 EEXS -= S*GeV; // New excitation energy (in MeV)
357
358 // NOTE: In-situ construction will be optimized away (no copying)
361
362 if (verboseLevel > 3)
363 G4cout << output.getOutgoingParticles().back() << G4endl;
364
365 ppout += mom;
366 } else {
367 if (itry_gam == itry_gam_max) try_again = false;
368 }
369 } // while (EEXS > cut_off
370 try_again = false;
371 } else { // if (prob_sum < prob_cut_off)
372 G4double SL = prob_sum * inuclRndm();
373 if (verboseLevel > 3) G4cout << " random SL " << SL << G4endl;
374
375 G4double S1 = 0.0;
376 for (i = 0; i < 7; i++) { // Select evaporation scenario
377 S1 += W[i];
378 if (SL <= S1) {
379 icase = i;
380 break;
381 };
382 };
383
384 if (icase < 0) continue; // Failed to choose scenario, try again
385
386 if (icase < 6) { // particle or light nuclei escape
387 if (verboseLevel > 2) {
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!")
392 << ")?" << G4endl;
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));
397 G4double d1 = 1.0 / ur;
398 G4double d2 = 1.0 / (ur - 1.0);
399 G4int itry1 = 0;
400 G4bool bad = true;
401
402 while (itry1 < itry_max && bad) {
403 itry1++;
404 G4int itry = 0;
405 G4double EPR = -1.0;
406 G4double S = 0.0;
407
408 while (itry < itry_max && EPR < 0.0) {
409 itry++;
410 G4double uu = uc + std::log((1.0 - d1) * inuclRndm() + d2);
411 S = 0.5 * (uc * uc - uu * uu) / u[icase];
412 EPR = TM[icase] - S * A / (A - 1.0) + V[icase];
413 };
414
415 if (verboseLevel > 3) {
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) { // real escape
421 if (verboseLevel > 2)
422 G4cout << " escape itry1 " << itry1 << " icase "
423 << icase << " S (MeV) " << S << G4endl;
424
425 S /= GeV; // Convert to Bertini units
426
427 if (icase < 2) { // particle escape
428 G4int ptype = 2 - icase;
429 if (verboseLevel > 2)
430 G4cout << " particle " << ptype << " escape" << G4endl;
431
432 // generate particle momentum
433 G4double mass =
435
436 G4double pmod = std::sqrt((2.0 * mass + S) * S);
438
439 // Push evaporated particle into current rest frame
440 mom = toTheNucleiSystemRestFrame.backToTheLab(mom);
441
442 if (verboseLevel > 2) {
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 // New excitation energy depends on residual nuclear state
450 G4double mass_new =
451 G4InuclNuclei::getNucleiMass(A1[icase],Z1[icase]);
452
453 G4double EEXS_new = ((PEX-mom).m() - mass_new)*GeV;
454 if (EEXS_new < 0.0) continue; // Sanity check for new nucleus
455
456 PEX -= mom; // Remaining four-momentum
457 EEXS = EEXS_new;
458
459 A = A1[icase];
460 Z = Z1[icase];
461
462 // NOTE: In-situ construction optimizes away (no copying)
465
466 if (verboseLevel > 3)
467 G4cout << output.getOutgoingParticles().back() << G4endl;
468
469 ppout += mom;
470 bad = false;
471 } else { // if (icase < 2)
472 if (verboseLevel > 2) {
473 G4cout << " nucleus A " << AN[icase] << " Z " << Q[icase]
474 << " escape icase " << icase << G4endl;
475 }
476
477 G4double mass =
478 G4InuclNuclei::getNucleiMass(AN[icase],Q[icase]);
479
480 // generate particle momentum
481 G4double pmod = std::sqrt((2.0 * mass + S) * S);
483
484 // Push evaporated particle into current rest frame
485 mom = toTheNucleiSystemRestFrame.backToTheLab(mom);
486
487 if (verboseLevel > 2) {
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 // New excitation energy depends on residual nuclear state
495 G4double mass_new =
496 G4InuclNuclei::getNucleiMass(A1[icase],Z1[icase]);
497
498 G4double EEXS_new = ((PEX-mom).m() - mass_new)*GeV;
499 if (EEXS_new < 0.0) continue; // Sanity check for new nucleus
500
501 PEX -= mom; // Remaining four-momentum
502 EEXS = EEXS_new;
503
504 A = A1[icase];
505 Z = Z1[icase];
506
507 // NOTE: In-situ constructor optimizes away (no copying)
509 AN[icase], Q[icase], 0.*GeV,
511
512 if (verboseLevel > 3)
513 G4cout << output.getOutgoingNuclei().back() << G4endl;
514
515 ppout += mom;
516 bad = false;
517 } // if (icase < 2)
518 } // if (EPR > 0.0 ...
519 } // while (itry1 ...
520
521 if (itry1 == itry_max || bad) try_again = false;
522 } else { // if (icase < 6)
524
525 if (verboseLevel > 2) {
526 G4cout << " fission: A " << A << " Z " << Z << " eexs " << EEXS
527 << " Wn " << W[0] << " Wf " << W[6] << G4endl;
528 }
529
530 // Catch fission output separately for verification
531 fission_output.reset();
532 theFissioner.collide(0, &nuclei, fission_output);
533
534 std::vector<G4InuclNuclei>& nuclea = fission_output.getOutgoingNuclei();
535 if (nuclea.size() == 2) { // fission ok
536 if (verboseLevel > 2) G4cout << " fission done in eql" << G4endl;
537
538 // Move fission fragments to lab frame for processing
539 fission_output.boostToLabFrame(toTheNucleiSystemRestFrame);
540
541 // Now evaporate the fission fragments individually
542 G4bool prevDoChecks = doConservationChecks; // Turn off checking
544
545 this->collide(0, &nuclea[0], output);
546 this->collide(0, &nuclea[1], output);
547
548 setConservationChecks(prevDoChecks); // Restore previous flag value
549 validateOutput(0, target, output); // Check energy conservation
550 return;
551 } else { // fission forbidden now
552 fission_open = false;
553 }
554 } // End of fission case
555 } // if (prob_sum < prob_cut_off)
556 } // while (try_again
557
558 // this time it's final nuclei
559
560 if (itry_global == itry_global_max) {
561 if (verboseLevel > 3) {
562 G4cout << " ! itry_global " << itry_global_max << G4endl;
563 }
564 }
565
566 G4LorentzVector pnuc = pin - ppout;
567
568 // NOTE: In-situ constructor will be optimized away (no copying)
569 output.addOutgoingNucleus(G4InuclNuclei(pnuc, A, Z, EEXS,
571
572 if (verboseLevel > 3) {
573 G4cout << " remaining nucleus \n" << output.getOutgoingNuclei().back()
574 << G4endl;
575 }
576
577
578 validateOutput(0, target, output); // Check energy conservation
579 return;
580}
@ photon
double G4double
Definition: G4Types.hh:64
int G4int
Definition: G4Types.hh:66
bool G4bool
Definition: G4Types.hh:67
#define G4endl
Definition: G4ios.hh:52
G4DLLIMPORT std::ostream G4cerr
G4DLLIMPORT std::ostream G4cout
void collide(G4InuclParticle *bullet, G4InuclParticle *target, G4CollisionOutput &output)
Definition: G4BigBanger.cc:65
virtual void setVerboseLevel(G4int verbose=0)
virtual void setConservationChecks(G4bool doBalance=true)
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)
Definition: G4Fissioner.cc:60
static G4double getParticleMass(G4int type)
G4double getNucleiMass() const
G4int getZ() const
G4double getExitationEnergy() const
G4int getA() const
G4LorentzVector getMomentum() const
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)
Definition: paraMaker.cc:38

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


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