66 using namespace G4InuclSpecialFunctions;
77 G4cout <<
" >>> G4NonEquilibriumEvaporator::deExcite" <<
G4endl;
82 const G4int a_cut = 5;
83 const G4int z_cut = 3;
88 const G4int itry_max = 1000;
101 G4int QP = QPP + QNP;
102 G4int QH = QPH + QNH;
108 toTheExitonSystemRestFrame.
setBullet(dummy);
117 G4bool try_again = (NEX > 0);
120 std::pair<G4double, G4double> parms;
123 if (
A >= a_cut &&
Z >= z_cut &&
EEXS > eexs_cut) {
131 G4cout <<
" A " <<
A <<
" Z " <<
Z <<
" mass " << nuc_mass
141 if (QEX < std::sqrt(2.0 * EG)) {
143 G4cout <<
" QEX " << QEX <<
" < sqrt(2*EG) " << std::sqrt(2.*EG)
144 <<
" NEX " << NEX <<
G4endl;
148 const G4double& CPA1 = parms.second;
160 G4cout <<
" AK1 " << AK1 <<
" CPA1 " <<
" VP " << VP
161 <<
"\n bind(A,Z) " << DM1 <<
" dBind(N) " << BN
162 <<
" dBind(P) " << BP
163 <<
"\n EMN " << EMN <<
" EMP " << EMP <<
G4endl;
166 if (EMN > eexs_cut) {
170 G4double APH = 0.25 * (QP * QP + QH * QH + QP - 3 * QH);
171 G4double APH1 = APH + 0.5 * (QP + QH);
176 G4cout <<
" APH " << APH <<
" APH1 " << APH1 <<
" ESP " << ESP
180 MELE *= std::sqrt(15.0 / ESP);
181 }
else if(ESP < 7.0) {
182 MELE *= std::sqrt(ESP / 7.0);
183 if (ESP < 2.0) MELE *= std::sqrt(ESP / 2.0);
190 G4cout <<
" MELE " << MELE <<
" F1 " << F1 <<
" F2 " << F2
193 if (F1 > 0.0 && F2 > 0.0) {
197 D[0] = M1 * F2 * F2 *
theG4Pow->
powN(F, NEX-1) / (QEX+1);
199 G4cout <<
" D[0] " << D[0] <<
" with F " << F
207 D[1] = 0.0462 / parlev /
G4cbrt(
A) * QP * EEXS / QEX;
210 D[2] = D[1] *
theG4Pow->
powN(EMP/EEXS, NEX) * (1.0 + CPA1);
214 G4cout <<
" D[1] " << D[1] <<
" with powN(EMN/EEXS, NEX) "
216 <<
" D[2] " << D[2] <<
" with powN(EMP/EEXS, NEX) "
220 if (QNP < 1) D[1] = 0.0;
221 if (QPP < 1) D[2] = 0.0;
223 try_again = NEX > 1 && (D[1] > width_cut * D[0] ||
224 D[2] > width_cut * D[0]);
234 for (
G4int i = 0; i < 3; i++) {
246 }
else try_again =
false;
247 }
else try_again =
false;
259 if (
A < 3.0) try_again =
false;
267 if (QNP < 1) icase = 0;
277 if (QPP < 1) icase = 0;
283 if (
Z-1 < 1) try_again =
false;
287 if (try_again && icase != 0) {
289 G4cout <<
" ptype " << ptype <<
" B " << B <<
" V " << V
295 if (E < 0.0) icase = 0;
304 while (itry1 < itry_max && icase > 0 && bad) {
309 while (EEXS_new < 0.0 && itry < itry_max) {
315 X = 1.0 - std::sqrt(R);
322 G4cout <<
" R " << R <<
" QEX2 " << QEX2
323 <<
" powA(R, QEX2) " << X <<
G4endl;
326 for (
G4int i = 0; i < 1000; i++) {
328 (1.0 + QEX2 * X * (1.0 - R /
theG4Pow->
powN(X, NEX)) / (1.0 - X));
330 G4cout <<
" NEX " << NEX <<
" powN(X, NEX) "
336 if (std::fabs(DX / X) < 0.01)
break;
341 EEXS_new = EB - EPART *
A / (
A-1);
344 if (itry == itry_max || EEXS_new < 0.0) {
350 G4cout <<
" particle " << ptype <<
" escape " <<
G4endl;
359 G4double pmod = std::sqrt(EPART * (2.0 * mass + EPART));
377 if (ptype == 2) QNP_new--;
381 <<
" pz " <<
PEX.
pz() <<
" E " <<
PEX.
e() << G4endl
382 <<
" evaporate px " << mom.
px() <<
" py " << mom.
py()
383 <<
" pz " << mom.
pz() <<
" E " << mom.
e() <<
G4endl;
389 EEXS_new = ((
PEX-mom).
m() - mass_new)*
GeV;
390 if (EEXS_new < 0.)
continue;
411 G4cout << particle << G4endl
412 <<
" ppout px " << ppout.
px() <<
" py " << ppout.
py()
413 <<
" pz " << ppout.
pz() <<
" E " << ppout.
e() <<
G4endl;
419 if (itry1 == itry_max) icase = 0;
425 if (icase == 0 && try_again) {
430 G4double XNUN = 1.0 / (1.6 + ESP / EFN);
431 G4double XNUP = 1.0 / (1.6 + ESP / EFP);
436 G4double PP = (QPP * SNN1 + QNP * SPN1) * ZR;
437 G4double PN = (QPP * SPN2 + QNP * SNN2) * (AR - ZR);
455 if (ZR < 2) try_again =
false;
457 }
else try_again =
false;
460 }
else try_again =
false;
461 }
else try_again =
false;
462 }
else try_again =
false;
483 G4cout <<
" >>> G4NonEquilibriumEvaporator::getMatrixElement" <<
G4endl;
488 if (a > 150) me = 100.0;
489 else if (a > 20) me = 140.0;
497 G4cout <<
" >>> G4NonEquilibriumEvaporator::getEO" <<
G4endl;
507 G4cout <<
" >>> G4NonEquilibriumEvaporator::getParLev" <<
G4endl;
void setModel(Model model)
std::vector< ExP01TrackerHit * > a
const G4Fragment & getRecoilFragment(G4int index=0) const
G4double getParLev(G4int A, G4int Z) const
void setBullet(const G4InuclParticle *bullet)
virtual G4bool validateOutput(const G4Fragment &target, G4CollisionOutput &output)
G4int protonQuasiParticles
G4LorentzVector generateWithRandomAngles(G4double p, G4double mass=0.)
void getTargetData(const G4Fragment &target)
G4double powN(G4double x, G4int n) const
void getTruncated(G4double Z, std::pair< G4double, G4double > &parms)
G4LorentzVector backToTheLab(const G4LorentzVector &mom) const
static constexpr double m
G4double csNN(G4double e)
void addRecoilFragment(const G4Fragment *aFragment)
G4double powA(G4double A, G4double y) const
G4NonEquilibriumEvaporator()
G4int neutronQuasiParticles
G4double G4cbrt(G4double x)
G4double getMatrixElement(G4int A) const
const G4Fragment & makeFragment(G4LorentzVector mom, G4int A, G4int Z, G4double EX=0.)
static const G4double * SL[nLA]
G4int numberOfOutgoingParticles() const
void addOutgoingParticle(const G4InuclElementaryParticle &particle)
void toTheTargetRestFrame()
std::map< G4String, G4AttDef > * GetInstance(const G4String &storeKey, G4bool &isNew)
void setTarget(const G4InuclParticle *target)
void setMomentum(const G4LorentzVector &mom)
G4double getNucleiMass() const
G4double getE0(G4int A) const
G4GLOB_DLL std::ostream G4cout
G4double csPN(G4double e)
virtual void deExcite(const G4Fragment &target, G4CollisionOutput &output)
G4double bindingEnergy(G4int A, G4int Z)
double B(double temperature)
static constexpr double GeV
G4double FermiEnergy(G4int A, G4int Z, G4int ntype)
void setVectM(const Hep3Vector &spatial, double mass)
G4InuclSpecialFunctions::paraMaker theParaMaker