G4EquilibriumEvaporator Class Reference

#include <G4EquilibriumEvaporator.hh>

Inheritance diagram for G4EquilibriumEvaporator:

G4CascadeColliderBase G4VCascadeCollider

Public Member Functions

 G4EquilibriumEvaporator ()
virtual ~G4EquilibriumEvaporator ()
void collide (G4InuclParticle *bullet, G4InuclParticle *target, G4CollisionOutput &output)

Detailed Description

Definition at line 49 of file G4EquilibriumEvaporator.hh.


Constructor & Destructor Documentation

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]

Definition at line 90 of file G4EquilibriumEvaporator.cc.

00090 {}


Member Function Documentation

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 }                    


The documentation for this class was generated from the following files:
Generated on Mon May 27 17:51:54 2013 for Geant4 by  doxygen 1.4.7