Geant4 9.6.0
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 ()
 
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
 
G4CascadeCheckBalance * balance
 
- Protected Attributes inherited from G4VCascadeCollider
const char * theName
 
G4int verboseLevel
 

Detailed Description

Definition at line 42 of file G4NonEquilibriumEvaporator.hh.

Constructor & Destructor Documentation

◆ G4NonEquilibriumEvaporator()

G4NonEquilibriumEvaporator::G4NonEquilibriumEvaporator ( )

Definition at line 62 of file G4NonEquilibriumEvaporator.cc.

63 : G4CascadeColliderBase("G4NonEquilibriumEvaporator") {}

◆ ~G4NonEquilibriumEvaporator()

virtual G4NonEquilibriumEvaporator::~G4NonEquilibriumEvaporator ( )
inlinevirtual

Definition at line 45 of file G4NonEquilibriumEvaporator.hh.

45{}

Member Function Documentation

◆ collide()

void G4NonEquilibriumEvaporator::collide ( G4InuclParticle *  bullet,
G4InuclParticle *  target,
G4CollisionOutput &  output 
)
virtual

Implements G4VCascadeCollider.

Definition at line 66 of file G4NonEquilibriumEvaporator.cc.

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

Referenced by G4CascadeDeexcitation::collide().

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