00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049 #include <cmath>
00050
00051 #include "G4NonEquilibriumEvaporator.hh"
00052 #include "G4SystemOfUnits.hh"
00053 #include "G4CollisionOutput.hh"
00054 #include "G4InuclElementaryParticle.hh"
00055 #include "G4InuclNuclei.hh"
00056 #include "G4InuclSpecialFunctions.hh"
00057 #include "G4LorentzConvertor.hh"
00058
00059 using namespace G4InuclSpecialFunctions;
00060
00061
00062 G4NonEquilibriumEvaporator::G4NonEquilibriumEvaporator()
00063 : G4CascadeColliderBase("G4NonEquilibriumEvaporator") {}
00064
00065
00066 void G4NonEquilibriumEvaporator::collide(G4InuclParticle* ,
00067 G4InuclParticle* target,
00068 G4CollisionOutput& output) {
00069
00070 if (verboseLevel) {
00071 G4cout << " >>> G4NonEquilibriumEvaporator::collide" << G4endl;
00072 }
00073
00074
00075 G4InuclNuclei* nuclei_target = dynamic_cast<G4InuclNuclei*>(target);
00076 if (!nuclei_target) {
00077 G4cerr << " NonEquilibriumEvaporator -> target is not nuclei " << G4endl;
00078 return;
00079 }
00080
00081 if (verboseLevel > 2) G4cout << " evaporating target:\n" << *target << G4endl;
00082
00083 const G4int a_cut = 5;
00084 const G4int z_cut = 3;
00085
00086 const G4double eexs_cut = 0.1;
00087
00088 const G4double coul_coeff = 1.4;
00089 const G4int itry_max = 1000;
00090 const G4double small_ekin = 1.0e-6;
00091 const G4double width_cut = 0.005;
00092
00093 G4int A = nuclei_target->getA();
00094 G4int Z = nuclei_target->getZ();
00095
00096 G4LorentzVector PEX = nuclei_target->getMomentum();
00097 G4LorentzVector pin = PEX;
00098
00099 G4double EEXS = nuclei_target->getExitationEnergy();
00100
00101 G4ExitonConfiguration config = nuclei_target->getExitonConfiguration();
00102
00103 G4int QPP = config.protonQuasiParticles;
00104 G4int QNP = config.neutronQuasiParticles;
00105 G4int QPH = config.protonHoles;
00106 G4int QNH = config.neutronHoles;
00107
00108 G4int QP = QPP + QNP;
00109 G4int QH = QPH + QNH;
00110 G4int QEX = QP + QH;
00111
00112 G4InuclElementaryParticle dummy(small_ekin, 1);
00113 G4LorentzConvertor toTheExitonSystemRestFrame;
00114
00115 toTheExitonSystemRestFrame.setBullet(dummy);
00116
00117 G4double EFN = FermiEnergy(A, Z, 0);
00118 G4double EFP = FermiEnergy(A, Z, 1);
00119
00120 G4int AR = A - QP;
00121 G4int ZR = Z - QPP;
00122 G4int NEX = QEX;
00123 G4LorentzVector ppout;
00124 G4bool try_again = (NEX > 0);
00125
00126
00127 std::pair<G4double, G4double> parms;
00128
00129 while (try_again) {
00130 if (A >= a_cut && Z >= z_cut && EEXS > eexs_cut) {
00131
00132 G4double nuc_mass = G4InuclNuclei::getNucleiMass(A, Z, EEXS);
00133 PEX.setVectM(PEX.vect(), nuc_mass);
00134 toTheExitonSystemRestFrame.setTarget(PEX);
00135 toTheExitonSystemRestFrame.toTheTargetRestFrame();
00136
00137 if (verboseLevel > 2) {
00138 G4cout << " A " << A << " Z " << Z << " mass " << nuc_mass
00139 << " EEXS " << EEXS << G4endl;
00140 }
00141
00142 G4double MEL = getMatrixElement(A);
00143 G4double E0 = getE0(A);
00144 G4double PL = getParLev(A, Z);
00145 G4double parlev = PL / A;
00146 G4double EG = PL * EEXS;
00147
00148 if (QEX < std::sqrt(2.0 * EG)) {
00149 if (verboseLevel > 3)
00150 G4cout << " QEX " << QEX << " < sqrt(2*EG) " << std::sqrt(2.*EG)
00151 << " NEX " << NEX << G4endl;
00152
00153 paraMakerTruncated(Z, parms);
00154 const G4double& AK1 = parms.first;
00155 const G4double& CPA1 = parms.second;
00156
00157 G4double VP = coul_coeff * Z * AK1 / (G4cbrt(A-1) + 1.0) /
00158 (1.0 + EEXS / E0);
00159 G4double DM1 = bindingEnergy(A,Z);
00160 G4double BN = DM1 - bindingEnergy(A-1,Z);
00161 G4double BP = DM1 - bindingEnergy(A-1,Z-1);
00162 G4double EMN = EEXS - BN;
00163 G4double EMP = EEXS - BP - VP * A / (A-1);
00164 G4double ESP = 0.0;
00165
00166 if (verboseLevel > 3) {
00167 G4cout << " AK1 " << AK1 << " CPA1 " << " VP " << VP
00168 << "\n bind(A,Z) " << DM1 << " dBind(N) " << BN
00169 << " dBind(P) " << BP
00170 << "\n EMN " << EMN << " EMP " << EMP << G4endl;
00171 }
00172
00173 if (EMN > eexs_cut) {
00174 G4int icase = 0;
00175
00176 if (NEX > 1) {
00177 G4double APH = 0.25 * (QP * QP + QH * QH + QP - 3 * QH);
00178 G4double APH1 = APH + 0.5 * (QP + QH);
00179 ESP = EEXS / QEX;
00180 G4double MELE = MEL / ESP / (A*A*A);
00181
00182 if (verboseLevel > 3)
00183 G4cout << " APH " << APH << " APH1 " << APH1 << " ESP " << ESP
00184 << G4endl;
00185
00186 if (ESP > 15.0) {
00187 MELE *= std::sqrt(15.0 / ESP);
00188 } else if(ESP < 7.0) {
00189 MELE *= std::sqrt(ESP / 7.0);
00190 if (ESP < 2.0) MELE *= std::sqrt(ESP / 2.0);
00191 };
00192
00193 G4double F1 = EG - APH;
00194 G4double F2 = EG - APH1;
00195
00196 if (verboseLevel > 3)
00197 G4cout << " MELE " << MELE << " F1 " << F1 << " F2 " << F2
00198 << G4endl;
00199
00200 if (F1 > 0.0 && F2 > 0.0) {
00201 G4double F = F2 / F1;
00202 G4double M1 = 2.77 * MELE * PL;
00203 G4double D[3] = { 0., 0., 0. };
00204 D[0] = M1 * F2 * F2 * std::pow(F, NEX-1) / (QEX+1);
00205
00206 if (D[0] > 0.0) {
00207
00208 if (NEX >= 2) {
00209 D[1] = 0.0462 / parlev / G4cbrt(A) * QP * EEXS / QEX;
00210
00211 if (EMP > eexs_cut)
00212 D[2] = D[1] * std::pow(EMP / EEXS, NEX) * (1.0 + CPA1);
00213 D[1] *= std::pow(EMN / EEXS, NEX) * getAL(A);
00214
00215 if (QNP < 1) D[1] = 0.0;
00216 if (QPP < 1) D[2] = 0.0;
00217
00218 try_again = NEX > 1 && (D[1] > width_cut * D[0] ||
00219 D[2] > width_cut * D[0]);
00220
00221 if (try_again) {
00222 G4double D5 = D[0] + D[1] + D[2];
00223 G4double SL = D5 * inuclRndm();
00224 G4double S1 = 0.;
00225
00226 if (verboseLevel > 3)
00227 G4cout << " D5 " << D5 << " SL " << SL << G4endl;
00228
00229 for (G4int i = 0; i < 3; i++) {
00230 S1 += D[i];
00231 if (SL <= S1) {
00232 icase = i;
00233 break;
00234 }
00235 }
00236
00237 if (verboseLevel > 3)
00238 G4cout << " got icase " << icase << G4endl;
00239 }
00240 }
00241 } else try_again = false;
00242 } else try_again = false;
00243 }
00244
00245 if (try_again) {
00246 if (icase > 0) {
00247 if (verboseLevel > 3)
00248 G4cout << " try_again icase " << icase << G4endl;
00249
00250 G4double V = 0.0;
00251 G4int ptype = 0;
00252 G4double B = 0.0;
00253
00254 if (A < 3.0) try_again = false;
00255
00256 if (try_again) {
00257
00258 if (icase == 1) {
00259 if (verboseLevel > 3)
00260 G4cout << " trying neutron escape" << G4endl;
00261
00262 if (QNP < 1) icase = 0;
00263 else {
00264 B = BN;
00265 V = 0.0;
00266 ptype = 2;
00267 };
00268 } else {
00269 if (verboseLevel > 3)
00270 G4cout << " trying proton escape" << G4endl;
00271
00272 if (QPP < 1) icase = 0;
00273 else {
00274 B = BP;
00275 V = VP;
00276 ptype = 1;
00277
00278 if (Z-1 < 1) try_again = false;
00279 };
00280 };
00281
00282 if (try_again && icase != 0) {
00283 if (verboseLevel > 3)
00284 G4cout << " ptype " << ptype << " B " << B << " V " << V
00285 << G4endl;
00286
00287 G4double EB = EEXS - B;
00288 G4double E = EB - V * A / (A-1);
00289
00290 if (E < 0.0) icase = 0;
00291 else {
00292 G4double E1 = EB - V;
00293 G4double EEXS_new = -1.;
00294 G4double EPART = 0.0;
00295 G4int itry1 = 0;
00296 G4bool bad = true;
00297
00298 while (itry1 < itry_max && icase > 0 && bad) {
00299 itry1++;
00300 G4int itry = 0;
00301
00302 while (EEXS_new < 0.0 && itry < itry_max) {
00303 itry++;
00304 G4double R = inuclRndm();
00305 G4double X;
00306
00307 if (NEX == 2) {
00308 X = 1.0 - std::sqrt(R);
00309
00310 } else {
00311 G4double QEX2 = 1.0 / QEX;
00312 G4double QEX1 = 1.0 / (QEX-1);
00313 X = std::pow(0.5 * R, QEX2);
00314
00315 for (G4int i = 0; i < 1000; i++) {
00316 G4double DX = X * QEX1 *
00317 (1.0 + QEX2 * X * (1.0 - R / std::pow(X, NEX)) / (1.0 - X));
00318 X -= DX;
00319
00320 if (std::fabs(DX / X) < 0.01) break;
00321
00322 };
00323 };
00324 EPART = EB - X * E1;
00325 EEXS_new = EB - EPART * A / (A-1);
00326 }
00327
00328 if (itry == itry_max || EEXS_new < 0.0) {
00329 icase = 0;
00330 continue;
00331 }
00332
00333 if (verboseLevel > 2)
00334 G4cout << " particle " << ptype << " escape " << G4endl;
00335
00336 EPART /= GeV;
00337
00338 G4InuclElementaryParticle particle(ptype);
00339 particle.setModel(G4InuclParticle::NonEquilib);
00340
00341
00342 G4double mass = particle.getMass();
00343 G4double pmod = std::sqrt(EPART * (2.0 * mass + EPART));
00344 G4LorentzVector mom = generateWithRandomAngles(pmod,mass);
00345
00346
00347 mom = toTheExitonSystemRestFrame.backToTheLab(mom);
00348
00349
00350 G4int QPP_new = QPP;
00351 G4int QNP_new = QNP;
00352
00353 G4int A_new = A-1;
00354 G4int Z_new = Z;
00355
00356 if (ptype == 1) {
00357 QPP_new--;
00358 Z_new--;
00359 };
00360
00361 if (ptype == 2) QNP_new--;
00362
00363 if (verboseLevel > 3) {
00364 G4cout << " nucleus px " << PEX.px() << " py " << PEX.py()
00365 << " pz " << PEX.pz() << " E " << PEX.e() << G4endl
00366 << " evaporate px " << mom.px() << " py " << mom.py()
00367 << " pz " << mom.pz() << " E " << mom.e() << G4endl;
00368 }
00369
00370
00371 G4double mass_new = G4InuclNuclei::getNucleiMass(A_new, Z_new);
00372
00373 EEXS_new = ((PEX-mom).m() - mass_new)*GeV;
00374 if (EEXS_new < 0.) continue;
00375
00376 if (verboseLevel > 3)
00377 G4cout << " EEXS_new " << EEXS_new << G4endl;
00378
00379 PEX -= mom;
00380 EEXS = EEXS_new;
00381
00382 A = A_new;
00383 Z = Z_new;
00384
00385 NEX--;
00386 QEX--;
00387 QP--;
00388 QPP = QPP_new;
00389 QNP = QNP_new;
00390
00391 particle.setMomentum(mom);
00392 output.addOutgoingParticle(particle);
00393 ppout += mom;
00394 if (verboseLevel > 3) {
00395 G4cout << particle << G4endl
00396 << " ppout px " << ppout.px() << " py " << ppout.py()
00397 << " pz " << ppout.pz() << " E " << ppout.e() << G4endl;
00398 }
00399
00400 bad = false;
00401 }
00402
00403 if (itry1 == itry_max) icase = 0;
00404 }
00405 }
00406 }
00407 }
00408
00409 if (icase == 0 && try_again) {
00410 if (verboseLevel > 3) G4cout << " adding excitons" << G4endl;
00411
00412 G4double TNN = 1.6 * EFN + ESP;
00413 G4double TNP = 1.6 * EFP + ESP;
00414 G4double XNUN = 1.0 / (1.6 + ESP / EFN);
00415 G4double XNUP = 1.0 / (1.6 + ESP / EFP);
00416 G4double SNN1 = csNN(TNP) * XNUP;
00417 G4double SNN2 = csNN(TNN) * XNUN;
00418 G4double SPN1 = csPN(TNP) * XNUP;
00419 G4double SPN2 = csPN(TNN) * XNUN;
00420 G4double PP = (QPP * SNN1 + QNP * SPN1) * ZR;
00421 G4double PN = (QPP * SPN2 + QNP * SNN2) * (AR - ZR);
00422 G4double PW = PP + PN;
00423 NEX += 2;
00424 QEX += 2;
00425 QP++;
00426 QH++;
00427 AR--;
00428
00429 if (AR > 1) {
00430 G4double SL = PW * inuclRndm();
00431
00432 if (SL > PP) {
00433 QNP++;
00434 QNH++;
00435 } else {
00436 QPP++;
00437 QPH++;
00438 ZR--;
00439 if (ZR < 2) try_again = false;
00440 }
00441 } else try_again = false;
00442 }
00443 }
00444 } else try_again = false;
00445 } else try_again = false;
00446 } else try_again = false;
00447 }
00448
00449
00450
00451 if (output.numberOfOutgoingParticles() == 0) {
00452 output.addOutgoingNucleus(*nuclei_target);
00453 } else {
00454 G4LorentzVector pnuc = pin - ppout;
00455 G4InuclNuclei nuclei(pnuc, A, Z, EEXS, G4InuclParticle::NonEquilib);
00456
00457 if (verboseLevel > 3) G4cout << " remaining nucleus\n" << nuclei << G4endl;
00458 output.addOutgoingNucleus(nuclei);
00459 }
00460
00461 validateOutput(0, target, output);
00462 return;
00463 }
00464
00465 G4double G4NonEquilibriumEvaporator::getMatrixElement(G4int A) const {
00466
00467 if (verboseLevel > 3) {
00468 G4cout << " >>> G4NonEquilibriumEvaporator::getMatrixElement" << G4endl;
00469 }
00470
00471 G4double me;
00472
00473 if (A > 150) me = 100.0;
00474 else if (A > 20) me = 140.0;
00475 else me = 70.0;
00476
00477 return me;
00478 }
00479
00480 G4double G4NonEquilibriumEvaporator::getE0(G4int ) const {
00481 if (verboseLevel > 3) {
00482 G4cout << " >>> G4NonEquilibriumEvaporator::getEO" << G4endl;
00483 }
00484
00485 const G4double e0 = 200.0;
00486
00487 return e0;
00488 }
00489
00490 G4double G4NonEquilibriumEvaporator::getParLev(G4int A,
00491 G4int ) const {
00492
00493 if (verboseLevel > 3) {
00494 G4cout << " >>> G4NonEquilibriumEvaporator::getParLev" << G4endl;
00495 }
00496
00497
00498 G4double pl = 0.125 * A;
00499
00500 return pl;
00501 }