65 using namespace G4InuclSpecialFunctions;
76 G4cout <<
" >>> G4NonEquilibriumEvaporator::deExcite" <<
G4endl;
81 const G4int a_cut = 5;
82 const G4int z_cut = 3;
87 const G4int itry_max = 1000;
100 G4int QP = QPP + QNP;
101 G4int QH = QPH + QNH;
107 toTheExitonSystemRestFrame.
setBullet(dummy);
116 G4bool try_again = (NEX > 0);
119 std::pair<G4double, G4double> parms;
122 if (
A >= a_cut &&
Z >= z_cut &&
EEXS > eexs_cut) {
130 G4cout <<
" A " <<
A <<
" Z " <<
Z <<
" mass " << nuc_mass
140 if (QEX < std::sqrt(2.0 * EG)) {
142 G4cout <<
" QEX " << QEX <<
" < sqrt(2*EG) " << std::sqrt(2.*EG)
143 <<
" NEX " << NEX <<
G4endl;
147 const G4double& CPA1 = parms.second;
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;
165 if (EMN > eexs_cut) {
169 G4double APH = 0.25 * (QP * QP + QH * QH + QP - 3 * QH);
170 G4double APH1 = APH + 0.5 * (QP + QH);
175 G4cout <<
" APH " << APH <<
" APH1 " << APH1 <<
" ESP " << ESP
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);
189 G4cout <<
" MELE " << MELE <<
" F1 " << F1 <<
" F2 " << F2
192 if (F1 > 0.0 && F2 > 0.0) {
196 D[0] = M1 * F2 * F2 * theG4Pow->
powN(F, NEX-1) / (QEX+1);
198 G4cout <<
" D[0] " << D[0] <<
" with F " << F
199 <<
" powN(F,NEX-1) " << theG4Pow->
powN(F, NEX-1)
206 D[1] = 0.0462 / parlev /
G4cbrt(
A) * QP * EEXS / QEX;
209 D[2] = D[1] * theG4Pow->
powN(EMP/EEXS, NEX) * (1.0 + CPA1);
210 D[1] *= theG4Pow->
powN(EMN/EEXS, NEX) *
getAL(
A);
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) "
219 if (QNP < 1) D[1] = 0.0;
220 if (QPP < 1) D[2] = 0.0;
222 try_again = NEX > 1 && (D[1] > width_cut * D[0] ||
223 D[2] > width_cut * D[0]);
233 for (
G4int i = 0; i < 3; i++) {
245 }
else try_again =
false;
246 }
else try_again =
false;
258 if (
A < 3.0) try_again =
false;
266 if (QNP < 1) icase = 0;
276 if (QPP < 1) icase = 0;
282 if (
Z-1 < 1) try_again =
false;
286 if (try_again && icase != 0) {
288 G4cout <<
" ptype " << ptype <<
" B " << B <<
" V " << V
294 if (E < 0.0) icase = 0;
302 while (itry1 < itry_max && icase > 0 && bad) {
306 while (EEXS_new < 0.0 && itry < itry_max) {
312 X = 1.0 - std::sqrt(R);
317 X = theG4Pow->
powA(0.5*R, QEX2);
319 G4cout <<
" R " << R <<
" QEX2 " << QEX2
320 <<
" powA(R, QEX2) " << X <<
G4endl;
323 for (
G4int i = 0; i < 1000; i++) {
325 (1.0 + QEX2 * X * (1.0 - R / theG4Pow->
powN(X, NEX)) / (1.0 - X));
327 G4cout <<
" NEX " << NEX <<
" powN(X, NEX) "
333 if (std::fabs(DX / X) < 0.01)
break;
338 EEXS_new = EB - EPART *
A / (
A-1);
341 if (itry == itry_max || EEXS_new < 0.0) {
347 G4cout <<
" particle " << ptype <<
" escape " <<
G4endl;
356 G4double pmod = std::sqrt(EPART * (2.0 * mass + EPART));
374 if (ptype == 2) QNP_new--;
378 <<
" pz " <<
PEX.
pz() <<
" E " <<
PEX.
e() << G4endl
379 <<
" evaporate px " << mom.
px() <<
" py " << mom.
py()
380 <<
" pz " << mom.
pz() <<
" E " << mom.
e() <<
G4endl;
386 EEXS_new = ((
PEX-mom).
m() - mass_new)*
GeV;
387 if (EEXS_new < 0.)
continue;
408 G4cout << particle << G4endl
409 <<
" ppout px " << ppout.
px() <<
" py " << ppout.
py()
410 <<
" pz " << ppout.
pz() <<
" E " << ppout.
e() <<
G4endl;
416 if (itry1 == itry_max) icase = 0;
422 if (icase == 0 && try_again) {
427 G4double XNUN = 1.0 / (1.6 + ESP / EFN);
428 G4double XNUP = 1.0 / (1.6 + ESP / EFP);
433 G4double PP = (QPP * SNN1 + QNP * SPN1) * ZR;
434 G4double PN = (QPP * SPN2 + QNP * SNN2) * (AR - ZR);
452 if (ZR < 2) try_again =
false;
454 }
else try_again =
false;
457 }
else try_again =
false;
458 }
else try_again =
false;
459 }
else try_again =
false;
478 G4double G4NonEquilibriumEvaporator::getMatrixElement(
G4int a)
const {
480 G4cout <<
" >>> G4NonEquilibriumEvaporator::getMatrixElement" <<
G4endl;
485 if (a > 150) me = 100.0;
486 else if (a > 20) me = 140.0;
494 G4cout <<
" >>> G4NonEquilibriumEvaporator::getEO" <<
G4endl;
504 G4cout <<
" >>> G4NonEquilibriumEvaporator::getParLev" <<
G4endl;
G4double powA(G4double A, G4double y) const
G4double powN(G4double x, G4int n) const
virtual void deExcite(const G4Fragment &target, G4CollisionOutput &output)
G4double csNN(G4double e)
G4NonEquilibriumEvaporator()
void addOutgoingParticle(const G4InuclElementaryParticle &particle)
void setBullet(const G4InuclParticle *bullet)
G4LorentzVector backToTheLab(const G4LorentzVector &mom) const
void setVectM(const Hep3Vector &spatial, double mass)
void getTruncated(G4double Z, std::pair< G4double, G4double > &parms)
void getTargetData(const G4Fragment &target)
G4GLOB_DLL std::ostream G4cout
G4double getNucleiMass() const
G4int numberOfOutgoingParticles() const
G4double csPN(G4double e)
G4LorentzVector generateWithRandomAngles(G4double p, G4double mass=0.)
G4double G4cbrt(G4double x)
const G4Fragment & makeFragment(G4LorentzVector mom, G4int A, G4int Z, G4double EX=0.)
void setModel(Model model)
virtual G4bool validateOutput(const G4Fragment &target, G4CollisionOutput &output)
void addRecoilFragment(const G4Fragment *aFragment)
const G4Fragment & getRecoilFragment(G4int index=0) const
void setMomentum(const G4LorentzVector &mom)
void toTheTargetRestFrame()
G4double bindingEnergy(G4int A, G4int Z)
G4double FermiEnergy(G4int A, G4int Z, G4int ntype)
std::map< G4String, G4AttDef > * GetInstance(G4String storeKey, G4bool &isNew)
G4int neutronQuasiParticles
void setTarget(const G4InuclParticle *target)
G4int protonQuasiParticles