Geant4 11.1.1
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
Public Member Functions | List of all members
G4NonEquilibriumEvaporator Class Reference

#include <G4NonEquilibriumEvaporator.hh>

+ Inheritance diagram for G4NonEquilibriumEvaporator:

Public Member Functions

 G4NonEquilibriumEvaporator ()
 
virtual ~G4NonEquilibriumEvaporator ()
 
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 G4Fragment & makeFragment (G4LorentzVector mom, G4int A, G4int Z, G4double EX=0.)
 
const G4Fragment & makeFragment (G4int A, G4int Z, G4double EX=0.)
 
- Protected Member Functions inherited from G4VCascadeCollider
virtual void setName (const G4String &name)
 
- Protected Attributes inherited from G4CascadeDeexciteBase
G4CascadeCheckBalance * balance
 
G4int A
 
G4int Z
 
G4LorentzVector PEX
 
G4double EEXS
 
G4Fragment aFragment
 
- Protected Attributes inherited from G4VCascadeCollider
G4String theName
 
G4int verboseLevel
 

Detailed Description

Definition at line 46 of file G4NonEquilibriumEvaporator.hh.

Constructor & Destructor Documentation

◆ G4NonEquilibriumEvaporator()

G4NonEquilibriumEvaporator::G4NonEquilibriumEvaporator ( )

Definition at line 68 of file G4NonEquilibriumEvaporator.cc.

69 : G4CascadeDeexciteBase("G4NonEquilibriumEvaporator"),
70 theParaMaker(verboseLevel), theG4Pow(G4Pow::GetInstance()) {}
G4Pow::GetInstance
static G4Pow * GetInstance()
Definition: G4Pow.cc:41

◆ ~G4NonEquilibriumEvaporator()

virtual G4NonEquilibriumEvaporator::~G4NonEquilibriumEvaporator ( )
inlinevirtual

Definition at line 49 of file G4NonEquilibriumEvaporator.hh.

49{}

Member Function Documentation

◆ deExcite()

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

Implements G4VCascadeDeexcitation.

Definition at line 73 of file G4NonEquilibriumEvaporator.cc.

74 {
75 if (verboseLevel) {
76 G4cout << " >>> G4NonEquilibriumEvaporator::deExcite" << G4endl;
77 }
78
79 if (verboseLevel>1) G4cout << " evaporating target:\n" << target << G4endl;
80
81 const G4int a_cut = 5;
82 const G4int z_cut = 3;
83
84 const G4double eexs_cut = 0.1;
85
86 const G4double coul_coeff = 1.4;
87 const G4int itry_max = 1000;
88 const G4double small_ekin = 1.0e-6;
89 const G4double width_cut = 0.005;
90
91 getTargetData(target);
92 G4LorentzVector pin = PEX; // Save original four-vector for later
93
94 G4ExitonConfiguration config(target);
95 G4int QPP = config.protonQuasiParticles;
96 G4int QNP = config.neutronQuasiParticles;
97 G4int QPH = config.protonHoles;
98 G4int QNH = config.neutronHoles;
99
100 G4int QP = QPP + QNP;
101 G4int QH = QPH + QNH;
102 G4int QEX = QP + QH;
103
104 G4InuclElementaryParticle dummy(small_ekin, 1);
105 G4LorentzConvertor toTheExitonSystemRestFrame;
106 //*** toTheExitonSystemRestFrame.setVerbose(verboseLevel);
107 toTheExitonSystemRestFrame.setBullet(dummy);
108
109 G4double EFN = FermiEnergy(A, Z, 0);
110 G4double EFP = FermiEnergy(A, Z, 1);
111
112 G4int AR = A - QP;
113 G4int ZR = Z - QPP;
114 G4int NEX = QEX;
115 G4LorentzVector ppout;
116 G4bool try_again = (NEX > 0);
117
118 // Buffer for parameter sets
119 std::pair<G4double, G4double> parms;
120
121 while (try_again) { /* Loop checking 08.06.2015 MHK */
122 if (A >= a_cut && Z >= z_cut && EEXS > eexs_cut) { // ok
123 // update exiton system (include excitation energy!)
124 G4double nuc_mass = G4InuclNuclei::getNucleiMass(A, Z, EEXS);
125 PEX.setVectM(PEX.vect(), nuc_mass);
126 toTheExitonSystemRestFrame.setTarget(PEX);
127 toTheExitonSystemRestFrame.toTheTargetRestFrame();
128
129 if (verboseLevel > 2) {
130 G4cout << " A " << A << " Z " << Z << " mass " << nuc_mass
131 << " EEXS " << EEXS << G4endl;
132 }
133
134 G4double MEL = getMatrixElement(A);
135 G4double E0 = getE0(A);
136 G4double PL = getParLev(A, Z);
137 G4double parlev = PL / A;
138 G4double EG = PL * EEXS;
139
140 if (QEX < std::sqrt(2.0 * EG)) { // ok
141 if (verboseLevel > 3)
142 G4cout << " QEX " << QEX << " < sqrt(2*EG) " << std::sqrt(2.*EG)
143 << " NEX " << NEX << G4endl;
144
145 theParaMaker.getTruncated(Z, parms);
146 const G4double& AK1 = parms.first;
147 const G4double& CPA1 = parms.second;
148
149 G4double VP = coul_coeff * Z * AK1 / (G4cbrt(A-1) + 1.0) /
150 (1.0 + EEXS / E0);
151 G4double DM1 = bindingEnergy(A,Z);
152 G4double BN = DM1 - bindingEnergy(A-1,Z);
153 G4double BP = DM1 - bindingEnergy(A-1,Z-1);
154 G4double EMN = EEXS - BN;
155 G4double EMP = EEXS - BP - VP * A / (A-1);
156 G4double ESP = 0.0;
157
158 if (verboseLevel > 3) {
159 G4cout << " AK1 " << AK1 << " CPA1 " << " VP " << VP
160 << "\n bind(A,Z) " << DM1 << " dBind(N) " << BN
161 << " dBind(P) " << BP
162 << "\n EMN " << EMN << " EMP " << EMP << G4endl;
163 }
164
165 if (EMN > eexs_cut) { // ok
166 G4int icase = 0;
167
168 if (NEX > 1) {
169 G4double APH = 0.25 * (QP * QP + QH * QH + QP - 3 * QH);
170 G4double APH1 = APH + 0.5 * (QP + QH);
171 ESP = EEXS / QEX;
172 G4double MELE = MEL / ESP / (A*A*A);
173
174 if (verboseLevel > 3)
175 G4cout << " APH " << APH << " APH1 " << APH1 << " ESP " << ESP
176 << G4endl;
177
178 if (ESP > 15.0) {
179 MELE *= std::sqrt(15.0 / ESP);
180 } else if(ESP < 7.0) {
181 MELE *= std::sqrt(ESP / 7.0);
182 if (ESP < 2.0) MELE *= std::sqrt(ESP / 2.0);
183 };
184
185 G4double F1 = EG - APH;
186 G4double F2 = EG - APH1;
187
188 if (verboseLevel > 3)
189 G4cout << " MELE " << MELE << " F1 " << F1 << " F2 " << F2
190 << G4endl;
191
192 if (F1 > 0.0 && F2 > 0.0) {
193 G4double F = F2 / F1;
194 G4double M1 = 2.77 * MELE * PL;
195 G4double D[3] = { 0., 0., 0. };
196 D[0] = M1 * F2 * F2 * theG4Pow->powN(F, NEX-1) / (QEX+1);
197 if (verboseLevel > 3) {
198 G4cout << " D[0] " << D[0] << " with F " << F
199 << " powN(F,NEX-1) " << theG4Pow->powN(F, NEX-1)
200 << G4endl;
201 }
202
203 if (D[0] > 0.0) {
204
205 if (NEX >= 2) {
206 D[1] = 0.0462 / parlev / G4cbrt(A) * QP * EEXS / QEX;
207
208 if (EMP > eexs_cut)
209 D[2] = D[1] * theG4Pow->powN(EMP/EEXS, NEX) * (1.0 + CPA1);
210 D[1] *= theG4Pow->powN(EMN/EEXS, NEX) * getAL(A);
211
212 if (verboseLevel > 3) {
213 G4cout << " D[1] " << D[1] << " with powN(EMN/EEXS, NEX) "
214 << theG4Pow->powN(EMN/EEXS, NEX) << G4endl
215 << " D[2] " << D[2] << " with powN(EMP/EEXS, NEX) "
216 << theG4Pow->powN(EMP/EEXS, NEX) << G4endl;
217 }
218
219 if (QNP < 1) D[1] = 0.0;
220 if (QPP < 1) D[2] = 0.0;
221
222 try_again = NEX > 1 && (D[1] > width_cut * D[0] ||
223 D[2] > width_cut * D[0]);
224
225 if (try_again) {
226 G4double D5 = D[0] + D[1] + D[2];
227 G4double SL = D5 * inuclRndm();
228 G4double S1 = 0.;
229
230 if (verboseLevel > 3)
231 G4cout << " D5 " << D5 << " SL " << SL << G4endl;
232
233 for (G4int i = 0; i < 3; i++) {
234 S1 += D[i];
235 if (SL <= S1) {
236 icase = i;
237 break;
238 }
239 }
240
241 if (verboseLevel > 3)
242 G4cout << " got icase " << icase << G4endl;
243 } // if (try_again)
244 } // if (NEX >= 2)
245 } else try_again = false; // if (D[0] > 0)
246 } else try_again = false; // if (F1>0 && F2>0)
247 } // if (NEX > 1)
248
249 if (try_again) {
250 if (icase > 0) { // N -> N-1 with particle escape
251 if (verboseLevel > 3)
252 G4cout << " try_again icase " << icase << G4endl;
253
254 G4double V = 0.0;
255 G4int ptype = 0;
256 G4double B = 0.0;
257
258 if (A < 3.0) try_again = false;
259
260 if (try_again) {
261
262 if (icase == 1) { // neutron escape
263 if (verboseLevel > 3)
264 G4cout << " trying neutron escape" << G4endl;
265
266 if (QNP < 1) icase = 0;
267 else {
268 B = BN;
269 V = 0.0;
270 ptype = 2;
271 };
272 } else { // proton esape
273 if (verboseLevel > 3)
274 G4cout << " trying proton escape" << G4endl;
275
276 if (QPP < 1) icase = 0;
277 else {
278 B = BP;
279 V = VP;
280 ptype = 1;
281
282 if (Z-1 < 1) try_again = false;
283 };
284 };
285
286 if (try_again && icase != 0) {
287 if (verboseLevel > 3)
288 G4cout << " ptype " << ptype << " B " << B << " V " << V
289 << G4endl;
290
291 G4double EB = EEXS - B;
292 G4double E = EB - V * A / (A-1);
293
294 if (E < 0.0) icase = 0;
295 else {
296 G4double E1 = EB - V;
297 G4double EEXS_new = -1.;
298 G4double EPART = 0.0;
299 G4int itry1 = 0;
300 G4bool bad = true;
301
302 /* Loop checking 08.06.2015 MHK */
303 while (itry1 < itry_max && icase > 0 && bad) {
304 itry1++;
305 G4int itry = 0;
306
307 /* Loop checking 08.06.2015 MHK */
308 while (EEXS_new < 0.0 && itry < itry_max) {
309 itry++;
310 G4double R = inuclRndm();
311 G4double X;
312
313 if (NEX == 2) {
314 X = 1.0 - std::sqrt(R);
315
316 } else {
317 G4double QEX2 = 1.0 / QEX;
318 G4double QEX1 = 1.0 / (QEX-1);
319 X = theG4Pow->powA(0.5*R, QEX2);
320 if (verboseLevel > 3) {
321 G4cout << " R " << R << " QEX2 " << QEX2
322 << " powA(R, QEX2) " << X << G4endl;
323 }
324
325 for (G4int i = 0; i < 1000; i++) {
326 G4double DX = X * QEX1 *
327 (1.0 + QEX2 * X * (1.0 - R / theG4Pow->powN(X, NEX)) / (1.0 - X));
328 if (verboseLevel > 3) {
329 G4cout << " NEX " << NEX << " powN(X, NEX) "
330 << theG4Pow->powN(X, NEX) << G4endl;
331 }
332
333 X -= DX;
334
335 if (std::fabs(DX / X) < 0.01) break;
336
337 };
338 };
339 EPART = EB - X * E1;
340 EEXS_new = EB - EPART * A / (A-1);
341 } // while (EEXS_new < 0.0...
342
343 if (itry == itry_max || EEXS_new < 0.0) {
344 icase = 0;
345 continue;
346 }
347
348 if (verboseLevel > 2)
349 G4cout << " particle " << ptype << " escape " << G4endl;
350
351 EPART /= GeV; // From MeV to GeV
352
353 G4InuclElementaryParticle particle(ptype);
354 particle.setModel(G4InuclParticle::NonEquilib);
355
356 // generate particle momentum
357 G4double mass = particle.getMass();
358 G4double pmod = std::sqrt(EPART * (2.0 * mass + EPART));
359 G4LorentzVector mom = generateWithRandomAngles(pmod,mass);
360
361 // Push evaporated paricle into current rest frame
362 mom = toTheExitonSystemRestFrame.backToTheLab(mom);
363
364 // Adjust quasiparticle and nucleon counts
365 G4int QPP_new = QPP;
366 G4int QNP_new = QNP;
367
368 G4int A_new = A-1;
369 G4int Z_new = Z;
370
371 if (ptype == 1) {
372 QPP_new--;
373 Z_new--;
374 };
375
376 if (ptype == 2) QNP_new--;
377
378 if (verboseLevel > 3) {
379 G4cout << " nucleus px " << PEX.px() << " py " << PEX.py()
380 << " pz " << PEX.pz() << " E " << PEX.e() << G4endl
381 << " evaporate px " << mom.px() << " py " << mom.py()
382 << " pz " << mom.pz() << " E " << mom.e() << G4endl;
383 }
384
385 // New excitation energy depends on residual nuclear state
386 G4double mass_new = G4InuclNuclei::getNucleiMass(A_new, Z_new);
387
388 EEXS_new = ((PEX-mom).m() - mass_new)*GeV;
389 if (EEXS_new < 0.) continue; // Sanity check for new nucleus
390
391 if (verboseLevel > 3)
392 G4cout << " EEXS_new " << EEXS_new << G4endl;
393
394 PEX -= mom;
395 EEXS = EEXS_new;
396
397 A = A_new;
398 Z = Z_new;
399
400 NEX--;
401 QEX--;
402 QP--;
403 QPP = QPP_new;
404 QNP = QNP_new;
405
406 particle.setMomentum(mom);
407 output.addOutgoingParticle(particle);
408 ppout += mom;
409 if (verboseLevel > 3) {
410 G4cout << particle << G4endl
411 << " ppout px " << ppout.px() << " py " << ppout.py()
412 << " pz " << ppout.pz() << " E " << ppout.e() << G4endl;
413 }
414
415 bad = false;
416 } // while (itry1<itry_max && icase>0
417
418 if (itry1 == itry_max) icase = 0;
419 } // if (E < 0.) [else]
420 } // if (try_again && icase != 0)
421 } // if (try_again)
422 } // if (icase > 0)
423
424 if (icase == 0 && try_again) { // N -> N + 2
425 if (verboseLevel > 3) G4cout << " adding excitons" << G4endl;
426
427 G4double TNN = 1.6 * EFN + ESP;
428 G4double TNP = 1.6 * EFP + ESP;
429 G4double XNUN = 1.0 / (1.6 + ESP / EFN);
430 G4double XNUP = 1.0 / (1.6 + ESP / EFP);
431 G4double SNN1 = csNN(TNP) * XNUP;
432 G4double SNN2 = csNN(TNN) * XNUN;
433 G4double SPN1 = csPN(TNP) * XNUP;
434 G4double SPN2 = csPN(TNN) * XNUN;
435 G4double PP = (QPP * SNN1 + QNP * SPN1) * ZR;
436 G4double PN = (QPP * SPN2 + QNP * SNN2) * (AR - ZR);
437 G4double PW = PP + PN;
438 NEX += 2;
439 QEX += 2;
440 QP++;
441 QH++;
442 AR--;
443
444 if (AR > 1) {
445 G4double SL = PW * inuclRndm();
446
447 if (SL > PP) {
448 QNP++;
449 QNH++;
450 } else {
451 QPP++;
452 QPH++;
453 ZR--;
454 if (ZR < 2) try_again = false;
455 }
456 } else try_again = false;
457 } // if (icase==0 && try_again)
458 } // if (try_again)
459 } else try_again = false; // if (EMN > eexs_cut)
460 } else try_again = false; // if (QEX < sqrg(2*EG)
461 } else try_again = false; // if (A > a_cut ...
462 } // while (try_again)
463
464 // everything finished, set output fragment
465
466 if (output.numberOfOutgoingParticles() == 0) {
467 output.addRecoilFragment(target);
468 } else {
469 G4LorentzVector pnuc = pin - ppout;
470 output.addRecoilFragment(makeFragment(pnuc, A, Z, EEXS));
471
472 if (verboseLevel>3)
473 G4cout << " remaining nucleus\n" << output.getRecoilFragment() << G4endl;
474 }
475
476 validateOutput(target, output); // Check energy conservation, etc.
477 return;
478}
B
G4double B(G4double temperature)
Definition: G4DNAElectronHoleRecombination.cc:68
D
G4double D(G4double temp)
Definition: G4DNAElectronHoleRecombination.cc:90
G4double
double G4double
Definition: G4Types.hh:83
G4bool
bool G4bool
Definition: G4Types.hh:86
G4int
int G4int
Definition: G4Types.hh:85
G4endl
#define G4endl
Definition: G4ios.hh:57
G4cout
G4GLOB_DLL std::ostream G4cout
CLHEP::HepLorentzVector::setVectM
void setVectM(const Hep3Vector &spatial, double mass)
CLHEP::HepLorentzVector::vect
Hep3Vector vect() const
G4CascadeDeexciteBase::getTargetData
void getTargetData(const G4Fragment &target)
Definition: G4CascadeDeexciteBase.cc:66
G4CascadeDeexciteBase::makeFragment
const G4Fragment & makeFragment(G4LorentzVector mom, G4int A, G4int Z, G4double EX=0.)
Definition: G4CascadeDeexciteBase.cc:86
G4CascadeDeexciteBase::validateOutput
virtual G4bool validateOutput(const G4Fragment &target, G4CollisionOutput &output)
Definition: G4CascadeDeexciteBase.cc:130
G4CollisionOutput::numberOfOutgoingParticles
G4int numberOfOutgoingParticles() const
Definition: G4CollisionOutput.hh:125
G4CollisionOutput::addRecoilFragment
void addRecoilFragment(const G4Fragment *aFragment)
Definition: G4CollisionOutput.hh:99
G4CollisionOutput::addOutgoingParticle
void addOutgoingParticle(const G4InuclElementaryParticle &particle)
Definition: G4CollisionOutput.hh:80
G4CollisionOutput::getRecoilFragment
const G4Fragment & getRecoilFragment(G4int index=0) const
Definition: G4CollisionOutput.cc:115
G4InuclNuclei::getNucleiMass
G4double getNucleiMass() const
Definition: G4InuclNuclei.hh:139
G4InuclSpecialFunctions::paraMaker::getTruncated
void getTruncated(G4double Z, std::pair< G4double, G4double > &parms)
Definition: paraMaker.cc:96
G4LorentzConvertor::setBullet
void setBullet(const G4InuclParticle *bullet)
Definition: G4LorentzConvertor.cc:66
G4LorentzConvertor::backToTheLab
G4LorentzVector backToTheLab(const G4LorentzVector &mom) const
Definition: G4LorentzConvertor.cc:134
G4LorentzConvertor::setTarget
void setTarget(const G4InuclParticle *target)
Definition: G4LorentzConvertor.cc:70
G4Pow::powN
G4double powN(G4double x, G4int n) const
Definition: G4Pow.cc:162
G4Pow::powA
G4double powA(G4double A, G4double y) const
Definition: G4Pow.hh:230
config
struct config_s config
G4InuclSpecialFunctions::bindingEnergy
G4double bindingEnergy(G4int A, G4int Z)
Definition: bindingEnergy.cc:36
G4InuclSpecialFunctions::generateWithRandomAngles
G4LorentzVector generateWithRandomAngles(G4double p, G4double mass=0.)
Definition: G4InuclSpecialFunctions.cc:177
G4InuclSpecialFunctions::FermiEnergy
G4double FermiEnergy(G4int A, G4int Z, G4int ntype)
Definition: G4InuclSpecialFunctions.cc:108
config_s
Definition: deflate.c:119

Referenced by G4CascadeDeexcitation::deExcite().

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