#include <G4EquilibriumEvaporator.hh>
Inheritance diagram for G4EquilibriumEvaporator:
Public Member Functions | |
G4EquilibriumEvaporator () | |
virtual | ~G4EquilibriumEvaporator () |
void | collide (G4InuclParticle *bullet, G4InuclParticle *target, G4CollisionOutput &output) |
Definition at line 49 of file G4EquilibriumEvaporator.hh.
G4EquilibriumEvaporator::G4EquilibriumEvaporator | ( | ) |
Definition at line 84 of file G4EquilibriumEvaporator.cc.
00085 : G4CascadeColliderBase("G4EquilibriumEvaporator") { 00086 parms.first.resize(6,0.); 00087 parms.second.resize(6,0.); 00088 }
G4EquilibriumEvaporator::~G4EquilibriumEvaporator | ( | ) | [virtual] |
void G4EquilibriumEvaporator::collide | ( | G4InuclParticle * | bullet, | |
G4InuclParticle * | target, | |||
G4CollisionOutput & | output | |||
) | [virtual] |
Implements G4VCascadeCollider.
Definition at line 93 of file G4EquilibriumEvaporator.cc.
References G4CollisionOutput::addOutgoingNucleus(), G4CollisionOutput::addOutgoingParticle(), G4LorentzConvertor::backToTheLab(), G4InuclSpecialFunctions::bindingEnergy(), G4CollisionOutput::boostToLabFrame(), G4Fissioner::collide(), G4BigBanger::collide(), G4CascadeColliderBase::doConservationChecks, G4InuclParticle::Equilib, G4InuclSpecialFunctions::G4cbrt(), G4cerr, G4cout, G4endl, G4InuclSpecialFunctions::generateWithRandomAngles(), G4InuclNuclei::getA(), G4InuclSpecialFunctions::getAL(), G4InuclNuclei::getExitationEnergy(), G4InuclParticle::getMomentum(), G4InuclNuclei::getNucleiMass(), G4CollisionOutput::getOutgoingNuclei(), G4CollisionOutput::getOutgoingParticles(), G4InuclElementaryParticle::getParticleMass(), G4InuclNuclei::getZ(), G4InuclSpecialFunctions::inuclRndm(), G4InuclParticleNames::nuclei, G4InuclSpecialFunctions::paraMaker(), photon, G4InuclParticleNames::proton, G4CollisionOutput::reset(), G4LorentzConvertor::setBullet(), G4CascadeColliderBase::setConservationChecks(), G4LorentzConvertor::setTarget(), G4CascadeColliderBase::setVerboseLevel(), G4LorentzConvertor::toTheTargetRestFrame(), G4CascadeColliderBase::validateOutput(), and G4VCascadeCollider::verboseLevel.
Referenced by G4EvaporationInuclCollider::collide(), and G4CascadeDeexcitation::collide().
00095 { 00096 if (verboseLevel) { 00097 G4cout << " >>> G4EquilibriumEvaporator::collide" << G4endl; 00098 } 00099 00100 // Sanity check 00101 G4InuclNuclei* nuclei_target = dynamic_cast<G4InuclNuclei*>(target); 00102 if (!nuclei_target) { 00103 G4cerr << " EquilibriumEvaporator -> target is not nuclei " << G4endl; 00104 return; 00105 } 00106 00107 if (verboseLevel>1) G4cout << " evaporating target: \n" << *target << G4endl; 00108 00109 theFissioner.setVerboseLevel(verboseLevel); 00110 theBigBanger.setVerboseLevel(verboseLevel); 00111 00112 // simple implementation of the equilibium evaporation a la Dostrowski 00113 const G4double huge_num = 50.0; 00114 const G4double small = -50.0; 00115 const G4double prob_cut_off = 1.0e-15; 00116 const G4double Q1[6] = { 0.0, 0.0, 2.23, 8.49, 7.72, 28.3 }; 00117 const G4int AN[6] = { 1, 1, 2, 3, 3, 4 }; 00118 const G4int Q[6] = { 0, 1, 1, 1, 2, 2 }; 00119 const G4double G[6] = { 2.0, 2.0, 6.0, 6.0, 6.0, 4.0 }; 00120 const G4double BE = 0.0063; 00121 const G4double fisssion_cut = 1000.0; 00122 const G4double cut_off_energy = 0.1; 00123 00124 const G4double BF = 0.0242; 00125 const G4int itry_max = 1000; 00126 const G4int itry_global_max = 1000; 00127 const G4double small_ekin = 1.0e-6; 00128 const G4int itry_gam_max = 100; 00129 00130 G4double W[8], u[6], V[6], TM[6]; 00131 G4int A1[6], Z1[6]; 00132 00133 G4int A = nuclei_target->getA(); 00134 G4int Z = nuclei_target->getZ(); 00135 G4LorentzVector PEX = nuclei_target->getMomentum(); 00136 G4double EEXS = nuclei_target->getExitationEnergy(); 00137 00138 if (verboseLevel > 3) G4cout << " after noeq: eexs " << EEXS << G4endl; 00139 00140 G4InuclElementaryParticle dummy(small_ekin, proton); 00141 G4LorentzConvertor toTheNucleiSystemRestFrame; 00142 //*** toTheNucleiSystemRestFrame.setVerbose(verboseLevel); 00143 toTheNucleiSystemRestFrame.setBullet(dummy); 00144 00145 G4LorentzVector ppout; 00146 00147 // See if fragment should just be dispersed 00148 if (explosion(A, Z, EEXS)) { 00149 if (verboseLevel > 1) G4cout << " big bang in eql start " << G4endl; 00150 theBigBanger.collide(0, target, output); 00151 00152 validateOutput(0, target, output); // Check energy conservation 00153 return; 00154 } 00155 00156 // If nucleus is in ground state, no evaporation 00157 if (EEXS < cut_off_energy) { 00158 if (verboseLevel > 1) G4cout << " no energy for evaporation" << G4endl; 00159 output.addOutgoingNucleus(*nuclei_target); 00160 00161 validateOutput(0, target, output); // Check energy conservation 00162 return; 00163 } 00164 00165 // Initialize evaporation attempts 00166 G4double coul_coeff = (A >= 100.0) ? 1.4 : 1.2; 00167 00168 G4LorentzVector pin = PEX; // Save original target for testing 00169 00170 G4bool try_again = true; 00171 G4bool fission_open = true; 00172 G4int itry_global = 0; 00173 00174 while (try_again && itry_global < itry_global_max) { 00175 itry_global++; 00176 00177 // Set rest frame of current (recoiling) nucleus 00178 toTheNucleiSystemRestFrame.setTarget(PEX); 00179 toTheNucleiSystemRestFrame.toTheTargetRestFrame(); 00180 00181 if (verboseLevel > 2) { 00182 G4double nuc_mass = G4InuclNuclei::getNucleiMass(A, Z, EEXS); 00183 G4cout << " A " << A << " Z " << Z << " mass " << nuc_mass 00184 << " EEXS " << EEXS << G4endl; 00185 } 00186 00187 if (explosion(A, Z, EEXS)) { // big bang 00188 if (verboseLevel > 2) 00189 G4cout << " big bang in eql step " << itry_global << G4endl; 00190 00191 G4InuclNuclei nuclei(PEX, A, Z, EEXS, G4InuclParticle::Equilib); 00192 theBigBanger.collide(0, &nuclei, output); 00193 00194 validateOutput(0, target, output); // Check energy conservation 00195 return; 00196 } 00197 00198 if (EEXS < cut_off_energy) { // Evaporation not possible 00199 if (verboseLevel > 2) 00200 G4cout << " no energy for evaporation in eql step " << itry_global 00201 << G4endl; 00202 00203 try_again = false; 00204 break; 00205 } 00206 00207 // Normal evaporation chain 00208 G4double E0 = getE0(A); 00209 G4double parlev = getPARLEVDEN(A, Z); 00210 G4double u1 = parlev * A; 00211 00212 paraMaker(Z, parms); 00213 const std::vector<G4double>& AK = parms.first; 00214 const std::vector<G4double>& CPA = parms.second; 00215 00216 G4double DM0 = bindingEnergy(A,Z); 00217 G4int i(0); 00218 00219 for (i = 0; i < 6; i++) { 00220 A1[i] = A - AN[i]; 00221 Z1[i] = Z - Q[i]; 00222 u[i] = parlev * A1[i]; 00223 V[i] = 0.0; 00224 TM[i] = -0.1; 00225 00226 if (goodRemnant(A1[i], Z1[i])) { 00227 G4double QB = DM0 - bindingEnergy(A1[i],Z1[i]) - Q1[i]; 00228 V[i] = coul_coeff * Z * Q[i] * AK[i] / (1.0 + EEXS / E0) / 00229 (G4cbrt(A1[i]) + G4cbrt(AN[i])); 00230 TM[i] = EEXS - QB - V[i] * A / A1[i]; 00231 if (verboseLevel > 3) { 00232 G4cout << " i=" << i << " : QB " << QB << " u " << u[i] 00233 << " V " << V[i] << " TM " << TM[i] << G4endl; 00234 } 00235 } 00236 } 00237 00238 G4double ue = 2.0 * std::sqrt(u1 * EEXS); 00239 G4double prob_sum = 0.0; 00240 00241 W[0] = 0.0; 00242 if (TM[0] > cut_off_energy) { 00243 G4double AL = getAL(A); 00244 W[0] = BE * G4cbrt(A1[0]*A1[0]) * G[0] * AL; 00245 G4double TM1 = 2.0 * std::sqrt(u[0] * TM[0]) - ue; 00246 00247 if (TM1 > huge_num) TM1 = huge_num; 00248 else if (TM1 < small) TM1 = small; 00249 00250 W[0] *= std::exp(TM1); 00251 prob_sum += W[0]; 00252 } 00253 00254 for (i = 1; i < 6; i++) { 00255 W[i] = 0.0; 00256 if (TM[i] > cut_off_energy) { 00257 W[i] = BE * G4cbrt(A1[i]*A1[i]) * G[i] * (1.0 + CPA[i]); 00258 G4double TM1 = 2.0 * std::sqrt(u[i] * TM[i]) - ue; 00259 00260 if (TM1 > huge_num) TM1 = huge_num; 00261 else if (TM1 < small) TM1 = small; 00262 00263 W[i] *= std::exp(TM1); 00264 prob_sum += W[i]; 00265 } 00266 } 00267 00268 // fisson part 00269 W[6] = 0.0; 00270 if (A >= 100.0 && fission_open) { 00271 G4double X2 = Z * Z / A; 00272 G4double X1 = 1.0 - 2.0 * Z / A; 00273 G4double X = 0.019316 * X2 / (1.0 - 1.79 * X1 * X1); 00274 G4double EF = EEXS - getQF(X, X2, A, Z, EEXS); 00275 00276 if (EF > 0.0) { 00277 G4double AF = u1 * getAF(X, A, Z, EEXS); 00278 G4double TM1 = 2.0 * std::sqrt(AF * EF) - ue; 00279 00280 if (TM1 > huge_num) TM1 = huge_num; 00281 else if (TM1 < small) TM1 = small; 00282 00283 W[6] = BF * std::exp(TM1); 00284 if (W[6] > fisssion_cut*W[0]) W[6] = fisssion_cut*W[0]; 00285 00286 prob_sum += W[6]; 00287 } 00288 } 00289 00290 // again time to decide what next 00291 if (verboseLevel > 2) { 00292 G4cout << " Evaporation probabilities: sum = " << prob_sum 00293 << "\n\t n " << W[0] << " p " << W[1] << " d " << W[2] 00294 << " He3 " << W[3] << "\n\t t " << W[4] << " alpha " << W[5] 00295 << " fission " << W[6] << G4endl; 00296 } 00297 00298 G4int icase = -1; 00299 00300 if (prob_sum < prob_cut_off) { // photon emission chain 00301 G4double UCR0 = 2.5 + 150.0 / A; 00302 G4double T00 = 1.0 / (std::sqrt(u1 / UCR0) - 1.25 / UCR0); 00303 00304 if (verboseLevel > 3) 00305 G4cout << " UCR0 " << UCR0 << " T00 " << T00 << G4endl; 00306 00307 G4int itry_gam = 0; 00308 while (EEXS > cut_off_energy && try_again) { 00309 itry_gam++; 00310 G4int itry = 0; 00311 G4double T04 = 4.0 * T00; 00312 G4double FMAX; 00313 00314 if (T04 < EEXS) { 00315 FMAX = (T04*T04*T04*T04) * std::exp((EEXS - T04) / T00); 00316 } else { 00317 FMAX = EEXS*EEXS*EEXS*EEXS; 00318 }; 00319 00320 if (verboseLevel > 3) 00321 G4cout << " T04 " << T04 << " FMAX (EEXS^4) " << FMAX << G4endl; 00322 00323 G4double S(0), X1(0); 00324 while (itry < itry_max) { 00325 itry++; 00326 S = EEXS * inuclRndm(); 00327 X1 = (S*S*S*S) * std::exp((EEXS - S) / T00); 00328 00329 if (X1 > FMAX * inuclRndm()) break; 00330 }; 00331 00332 if (itry == itry_max) { // Maximum attempts exceeded 00333 try_again = false; 00334 break; 00335 } 00336 00337 if (verboseLevel > 2) G4cout << " photon escape ?" << G4endl; 00338 00339 if (S < EEXS) { // Valid evaporate 00340 S /= GeV; // Convert to Bertini units 00341 00342 G4double pmod = S; 00343 G4LorentzVector mom = generateWithRandomAngles(pmod, 0.); 00344 00345 // Push evaporated particle into current rest frame 00346 mom = toTheNucleiSystemRestFrame.backToTheLab(mom); 00347 00348 if (verboseLevel > 3) { 00349 G4cout << " nucleus px " << PEX.px() << " py " << PEX.py() 00350 << " pz " << PEX.pz() << " E " << PEX.e() << G4endl 00351 << " evaporate px " << mom.px() << " py " << mom.py() 00352 << " pz " << mom.pz() << " E " << mom.e() << G4endl; 00353 } 00354 00355 PEX -= mom; // Remaining four-momentum 00356 EEXS -= S*GeV; // New excitation energy (in MeV) 00357 00358 // NOTE: In-situ construction will be optimized away (no copying) 00359 output.addOutgoingParticle(G4InuclElementaryParticle(mom, photon, 00360 G4InuclParticle::Equilib)); 00361 00362 if (verboseLevel > 3) 00363 G4cout << output.getOutgoingParticles().back() << G4endl; 00364 00365 ppout += mom; 00366 } else { 00367 if (itry_gam == itry_gam_max) try_again = false; 00368 } 00369 } // while (EEXS > cut_off 00370 try_again = false; 00371 } else { // if (prob_sum < prob_cut_off) 00372 G4double SL = prob_sum * inuclRndm(); 00373 if (verboseLevel > 3) G4cout << " random SL " << SL << G4endl; 00374 00375 G4double S1 = 0.0; 00376 for (i = 0; i < 7; i++) { // Select evaporation scenario 00377 S1 += W[i]; 00378 if (SL <= S1) { 00379 icase = i; 00380 break; 00381 }; 00382 }; 00383 00384 if (icase < 0) continue; // Failed to choose scenario, try again 00385 00386 if (icase < 6) { // particle or light nuclei escape 00387 if (verboseLevel > 2) { 00388 G4cout << " particle/light-ion escape (" 00389 << (icase==0 ? "neutron" : icase==1 ? "proton" : 00390 icase==2 ? "deuteron" : icase==3 ? "He3" : 00391 icase==4 ? "triton" : icase==5 ? "alpha" : "ERROR!") 00392 << ")?" << G4endl; 00393 } 00394 00395 G4double uc = 2.0 * std::sqrt(u[icase] * TM[icase]); 00396 G4double ur = (uc > huge_num ? std::exp(huge_num) : std::exp(uc)); 00397 G4double d1 = 1.0 / ur; 00398 G4double d2 = 1.0 / (ur - 1.0); 00399 G4int itry1 = 0; 00400 G4bool bad = true; 00401 00402 while (itry1 < itry_max && bad) { 00403 itry1++; 00404 G4int itry = 0; 00405 G4double EPR = -1.0; 00406 G4double S = 0.0; 00407 00408 while (itry < itry_max && EPR < 0.0) { 00409 itry++; 00410 G4double uu = uc + std::log((1.0 - d1) * inuclRndm() + d2); 00411 S = 0.5 * (uc * uc - uu * uu) / u[icase]; 00412 EPR = TM[icase] - S * A / (A - 1.0) + V[icase]; 00413 }; 00414 00415 if (verboseLevel > 3) { 00416 G4cout << " EPR " << EPR << " V " << V[icase] 00417 << " S " << S << " EEXS " << EEXS << G4endl; 00418 } 00419 00420 if (EPR > 0.0 && S > V[icase] && S < EEXS) { // real escape 00421 if (verboseLevel > 2) 00422 G4cout << " escape itry1 " << itry1 << " icase " 00423 << icase << " S (MeV) " << S << G4endl; 00424 00425 S /= GeV; // Convert to Bertini units 00426 00427 if (icase < 2) { // particle escape 00428 G4int ptype = 2 - icase; 00429 if (verboseLevel > 2) 00430 G4cout << " particle " << ptype << " escape" << G4endl; 00431 00432 // generate particle momentum 00433 G4double mass = 00434 G4InuclElementaryParticle::getParticleMass(ptype); 00435 00436 G4double pmod = std::sqrt((2.0 * mass + S) * S); 00437 G4LorentzVector mom = generateWithRandomAngles(pmod, mass); 00438 00439 // Push evaporated particle into current rest frame 00440 mom = toTheNucleiSystemRestFrame.backToTheLab(mom); 00441 00442 if (verboseLevel > 2) { 00443 G4cout << " nucleus px " << PEX.px() << " py " << PEX.py() 00444 << " pz " << PEX.pz() << " E " << PEX.e() << G4endl 00445 << " evaporate px " << mom.px() << " py " << mom.py() 00446 << " pz " << mom.pz() << " E " << mom.e() << G4endl; 00447 } 00448 00449 // New excitation energy depends on residual nuclear state 00450 G4double mass_new = 00451 G4InuclNuclei::getNucleiMass(A1[icase],Z1[icase]); 00452 00453 G4double EEXS_new = ((PEX-mom).m() - mass_new)*GeV; 00454 if (EEXS_new < 0.0) continue; // Sanity check for new nucleus 00455 00456 PEX -= mom; // Remaining four-momentum 00457 EEXS = EEXS_new; 00458 00459 A = A1[icase]; 00460 Z = Z1[icase]; 00461 00462 // NOTE: In-situ construction optimizes away (no copying) 00463 output.addOutgoingParticle(G4InuclElementaryParticle(mom, 00464 ptype, G4InuclParticle::Equilib)); 00465 00466 if (verboseLevel > 3) 00467 G4cout << output.getOutgoingParticles().back() << G4endl; 00468 00469 ppout += mom; 00470 bad = false; 00471 } else { // if (icase < 2) 00472 if (verboseLevel > 2) { 00473 G4cout << " nucleus A " << AN[icase] << " Z " << Q[icase] 00474 << " escape icase " << icase << G4endl; 00475 } 00476 00477 G4double mass = 00478 G4InuclNuclei::getNucleiMass(AN[icase],Q[icase]); 00479 00480 // generate particle momentum 00481 G4double pmod = std::sqrt((2.0 * mass + S) * S); 00482 G4LorentzVector mom = generateWithRandomAngles(pmod,mass); 00483 00484 // Push evaporated particle into current rest frame 00485 mom = toTheNucleiSystemRestFrame.backToTheLab(mom); 00486 00487 if (verboseLevel > 2) { 00488 G4cout << " nucleus px " << PEX.px() << " py " << PEX.py() 00489 << " pz " << PEX.pz() << " E " << PEX.e() << G4endl 00490 << " evaporate px " << mom.px() << " py " << mom.py() 00491 << " pz " << mom.pz() << " E " << mom.e() << G4endl; 00492 } 00493 00494 // New excitation energy depends on residual nuclear state 00495 G4double mass_new = 00496 G4InuclNuclei::getNucleiMass(A1[icase],Z1[icase]); 00497 00498 G4double EEXS_new = ((PEX-mom).m() - mass_new)*GeV; 00499 if (EEXS_new < 0.0) continue; // Sanity check for new nucleus 00500 00501 PEX -= mom; // Remaining four-momentum 00502 EEXS = EEXS_new; 00503 00504 A = A1[icase]; 00505 Z = Z1[icase]; 00506 00507 // NOTE: In-situ constructor optimizes away (no copying) 00508 output.addOutgoingNucleus(G4InuclNuclei(mom, 00509 AN[icase], Q[icase], 0.*GeV, 00510 G4InuclParticle::Equilib)); 00511 00512 if (verboseLevel > 3) 00513 G4cout << output.getOutgoingNuclei().back() << G4endl; 00514 00515 ppout += mom; 00516 bad = false; 00517 } // if (icase < 2) 00518 } // if (EPR > 0.0 ... 00519 } // while (itry1 ... 00520 00521 if (itry1 == itry_max || bad) try_again = false; 00522 } else { // if (icase < 6) 00523 G4InuclNuclei nuclei(A, Z, EEXS, G4InuclParticle::Equilib); 00524 00525 if (verboseLevel > 2) { 00526 G4cout << " fission: A " << A << " Z " << Z << " eexs " << EEXS 00527 << " Wn " << W[0] << " Wf " << W[6] << G4endl; 00528 } 00529 00530 // Catch fission output separately for verification 00531 fission_output.reset(); 00532 theFissioner.collide(0, &nuclei, fission_output); 00533 00534 std::vector<G4InuclNuclei>& nuclea = fission_output.getOutgoingNuclei(); 00535 if (nuclea.size() == 2) { // fission ok 00536 if (verboseLevel > 2) G4cout << " fission done in eql" << G4endl; 00537 00538 // Move fission fragments to lab frame for processing 00539 fission_output.boostToLabFrame(toTheNucleiSystemRestFrame); 00540 00541 // Now evaporate the fission fragments individually 00542 G4bool prevDoChecks = doConservationChecks; // Turn off checking 00543 setConservationChecks(false); 00544 00545 this->collide(0, &nuclea[0], output); 00546 this->collide(0, &nuclea[1], output); 00547 00548 setConservationChecks(prevDoChecks); // Restore previous flag value 00549 validateOutput(0, target, output); // Check energy conservation 00550 return; 00551 } else { // fission forbidden now 00552 fission_open = false; 00553 } 00554 } // End of fission case 00555 } // if (prob_sum < prob_cut_off) 00556 } // while (try_again 00557 00558 // this time it's final nuclei 00559 00560 if (itry_global == itry_global_max) { 00561 if (verboseLevel > 3) { 00562 G4cout << " ! itry_global " << itry_global_max << G4endl; 00563 } 00564 } 00565 00566 G4LorentzVector pnuc = pin - ppout; 00567 00568 // NOTE: In-situ constructor will be optimized away (no copying) 00569 output.addOutgoingNucleus(G4InuclNuclei(pnuc, A, Z, EEXS, 00570 G4InuclParticle::Equilib)); 00571 00572 if (verboseLevel > 3) { 00573 G4cout << " remaining nucleus \n" << output.getOutgoingNuclei().back() 00574 << G4endl; 00575 } 00576 00577 00578 validateOutput(0, target, output); // Check energy conservation 00579 return; 00580 }