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;
   150     G4double VP = coul_coeff * 
Z * AK1 / (G4cbrt(
A-1) + 1.0) /
   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;  
   407               particle.setMomentum(mom);
   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;       
 
G4double powN(G4double x, G4int n) const
 
void addOutgoingParticle(const G4InuclElementaryParticle &particle)
 
void setBullet(const G4InuclParticle *bullet)
 
void setVectM(const Hep3Vector &spatial, double mass)
 
void getTargetData(const G4Fragment &target)
 
G4InuclSpecialFunctions::paraMaker theParaMaker
 
G4GLOB_DLL std::ostream G4cout
 
G4double getMatrixElement(G4int A) const
 
G4LorentzVector backToTheLab(const G4LorentzVector &mom) const
 
G4double getE0(G4int A) const
 
G4double getNucleiMass() const
 
G4double getParLev(G4int A, G4int Z) const
 
G4int numberOfOutgoingParticles() const
 
const G4Fragment & makeFragment(G4LorentzVector mom, G4int A, G4int Z, G4double EX=0.)
 
static const G4double * SL[nLA]
 
virtual G4bool validateOutput(const G4Fragment &target, G4CollisionOutput &output)
 
void addRecoilFragment(const G4Fragment *aFragment)
 
G4double powA(G4double A, G4double y) const
 
void toTheTargetRestFrame()
 
const G4Fragment & getRecoilFragment(G4int index=0) const
 
void setTarget(const G4InuclParticle *target)