G4EquilibriumEvaporator Class Reference

#include <G4EquilibriumEvaporator.hh>

Inheritance diagram for G4EquilibriumEvaporator:

Public Member Functions
	G4EquilibriumEvaporator ()
virtual	~G4EquilibriumEvaporator ()
virtual void	setVerboseLevel (G4int verbose)
virtual void	deExcite (const G4Fragment &target, G4CollisionOutput &output)
Public Member Functions inherited from G4CascadeDeexciteBase
	G4CascadeDeexciteBase (const char *name)
virtual	~G4CascadeDeexciteBase ()
Public Member Functions inherited from G4VCascadeDeexcitation
	G4VCascadeDeexcitation (const char *name)
virtual	~G4VCascadeDeexcitation ()
virtual void	collide (G4InuclParticle *bullet, G4InuclParticle *target, G4CollisionOutput &globalOutput)
Public Member Functions inherited from G4VCascadeCollider
	G4VCascadeCollider (const char *name, G4int verbose=0)
virtual	~G4VCascadeCollider ()

Additional Inherited Members
Protected Member Functions inherited from G4CascadeDeexciteBase
virtual G4bool	validateOutput (const G4Fragment &target, G4CollisionOutput &output)
virtual G4bool	validateOutput (const G4Fragment &target, const std::vector< G4InuclElementaryParticle > &particles)
virtual G4bool	validateOutput (const G4Fragment &target, const std::vector< G4InuclNuclei > &fragments)
void	getTargetData (const G4Fragment &target)
const G4Fragment &	makeFragment (G4LorentzVector mom, G4int A, G4int Z, G4double EX=0.)
const G4Fragment &	makeFragment (G4int A, G4int Z, G4double EX=0.)
Protected Member Functions inherited from G4VCascadeCollider
virtual void	setName (const char *name)
Protected Attributes inherited from G4CascadeDeexciteBase
G4CascadeCheckBalance *	balance
G4int	A
G4int	Z
G4LorentzVector	PEX
G4double	EEXS
G4Fragment	aFragment
Protected Attributes inherited from G4VCascadeCollider
const char *	theName
G4int	verboseLevel

Detailed Description

Definition at line 58 of file G4EquilibriumEvaporator.hh.

Constructor & Destructor Documentation

G4EquilibriumEvaporator::G4EquilibriumEvaporator

(

)

Definition at line 150 of file G4EquilibriumEvaporator.cc.

  : G4CascadeDeexciteBase("G4EquilibriumEvaporator"),
    theParaMaker(verboseLevel), QFinterp(XREP) {
  parms.first.resize(6,0.);
  parms.second.resize(6,0.);
}

G4EquilibriumEvaporator::~G4EquilibriumEvaporator ( )

virtual

Definition at line 157 of file G4EquilibriumEvaporator.cc.

{}

Member Function Documentation

void G4EquilibriumEvaporator::deExcite ( const G4Fragment &  target, G4CollisionOutput &  output  )

virtual

Implements G4VCascadeDeexcitation.

Definition at line 168 of file G4EquilibriumEvaporator.cc.

References G4CascadeDeexciteBase::A, G4CollisionOutput::addOutgoingNucleus(), G4CollisionOutput::addOutgoingParticle(), G4CollisionOutput::addRecoilFragment(), G4LorentzConvertor::backToTheLab(), G4InuclSpecialFunctions::bindingEnergy(), G4CollisionOutput::boostToLabFrame(), G4BigBanger::deExcite(), G4Fissioner::deExcite(), CLHEP::HepLorentzVector::e(), G4CascadeDeexciteBase::EEXS, G4InuclParticle::Equilib, G4InuclSpecialFunctions::G4cbrt(), G4cout, G4endl, G4Exp(), G4Log(), G4InuclSpecialFunctions::generateWithRandomAngles(), G4InuclSpecialFunctions::getAL(), G4InuclNuclei::getNucleiMass(), G4CollisionOutput::getOutgoingNuclei(), G4CollisionOutput::getOutgoingParticles(), G4InuclSpecialFunctions::paraMaker::getParams(), G4InuclElementaryParticle::getParticleMass(), G4CollisionOutput::getRecoilFragment(), G4CascadeDeexciteBase::getTargetData(), python.hepunit::GeV, G4InuclSpecialFunctions::inuclRndm(), python.hepunit::m, G4CascadeDeexciteBase::makeFragment(), G4CollisionOutput::numberOfFragments(), G4CascadeDeexciteBase::PEX, photon, G4InuclParticleNames::proton, CLHEP::HepLorentzVector::px(), CLHEP::HepLorentzVector::py(), CLHEP::HepLorentzVector::pz(), G4CollisionOutput::reset(), G4LorentzConvertor::setBullet(), G4LorentzConvertor::setTarget(), G4LorentzConvertor::toTheTargetRestFrame(), G4CascadeDeexciteBase::validateOutput(), G4VCascadeCollider::verboseLevel, and G4CascadeDeexciteBase::Z.

Referenced by G4EvaporationInuclCollider::deExcite(), and G4CascadeDeexcitation::deExcite().

                                                  {
  if (verboseLevel) {
    G4cout << " >>> G4EquilibriumEvaporator::deExcite" << G4endl;
  }

  if (verboseLevel>1) G4cout << " evaporating target: \n" << target << G4endl;

  // simple implementation of the equilibium evaporation a la Dostrowski
  const G4double huge_num = 50.0;
  const G4double small = -50.0;
  const G4double prob_cut_off = 1.0e-15;
  const G4double Q1[6] = { 0.0, 0.0, 2.23, 8.49, 7.72, 28.3 };
  const G4int AN[6] = { 1, 1, 2, 3, 3, 4 };
  const G4int Q[6] =  { 0, 1, 1, 1, 2, 2 };
  const G4double G[6] = { 2.0, 2.0, 6.0, 6.0, 6.0, 4.0 };
  const G4double BE = 0.0063;
  const G4double fisssion_cut = 1000.0;
  const G4double cut_off_energy = 0.1;

  const G4double BF = 0.0242;
  const G4int itry_max = 1000;
  const G4int itry_global_max = 1000;
  const G4double small_ekin = 1.0e-6;
  const G4int itry_gam_max = 100;

  G4double W[8], u[6], V[6], TM[6];
  G4int A1[6], Z1[6];

  getTargetData(target);
  if (verboseLevel > 3) G4cout << " after noeq: eexs " << EEXS << G4endl;

  G4InuclElementaryParticle dummy(small_ekin, proton);
  G4LorentzConvertor toTheNucleiSystemRestFrame;
  //*** toTheNucleiSystemRestFrame.setVerbose(verboseLevel);
  toTheNucleiSystemRestFrame.setBullet(dummy);

  G4LorentzVector ppout;
  
  // See if fragment should just be dispersed
  if (explosion(A, Z, EEXS)) {
    if (verboseLevel > 1) G4cout << " big bang in eql start " << G4endl;
    theBigBanger.deExcite(target, output);

    validateOutput(target, output);     // Check energy conservation
    return;
  }

  // If nucleus is in ground state, no evaporation
  if (EEXS < cut_off_energy) {
    if (verboseLevel > 1) G4cout << " no energy for evaporation" << G4endl;
    output.addRecoilFragment(target);

    validateOutput(target, output);     // Check energy conservation
    return;
  }

  // Initialize evaporation attempts
  G4double coul_coeff = (A >= 100.0) ? 1.4 : 1.2;
   
  G4LorentzVector pin = PEX;    // Save original target for testing
    
  G4bool try_again = true;  
  G4bool fission_open = true;
  G4int itry_global = 0;
    
  while (try_again && itry_global < itry_global_max) {
    itry_global++;

    // Set rest frame of current (recoiling) nucleus
    toTheNucleiSystemRestFrame.setTarget(PEX);
    toTheNucleiSystemRestFrame.toTheTargetRestFrame();
      
    if (verboseLevel > 2) {
      G4double nuc_mass = G4InuclNuclei::getNucleiMass(A, Z, EEXS);
      G4cout << " A " << A << " Z " << Z << " mass " << nuc_mass
         << " EEXS " << EEXS << G4endl;
    }
      
    if (explosion(A, Z, EEXS)) {            // big bang
      if (verboseLevel > 2) 
    G4cout << " big bang in eql step " << itry_global << G4endl;

      theBigBanger.deExcite(makeFragment(PEX,A,Z,EEXS), output);

      validateOutput(target, output);   // Check energy conservation
      return;   
    } 

    if (EEXS < cut_off_energy) {    // Evaporation not possible
      if (verboseLevel > 2)
    G4cout << " no energy for evaporation in eql step " << itry_global 
           << G4endl;

      try_again = false;
      break;
    }

    // Normal evaporation chain
    G4double E0 = getE0(A);
    G4double parlev = getPARLEVDEN(A, Z);
    G4double u1 = parlev * A;

    theParaMaker.getParams(Z, parms);
    const std::vector<G4double>& AK = parms.first;
    const std::vector<G4double>& CPA = parms.second;

    G4double DM0 = bindingEnergy(A,Z);
    G4int i(0);

    for (i = 0; i < 6; i++) {
      A1[i] = A - AN[i];
      Z1[i] = Z - Q[i];
      u[i] = parlev * A1[i];
      V[i] = 0.0;
      TM[i] = -0.1;

      if (goodRemnant(A1[i], Z1[i])) {
    G4double QB = DM0 - bindingEnergy(A1[i],Z1[i]) - Q1[i];
    V[i] = coul_coeff * Z * Q[i] * AK[i] / (1.0 + EEXS / E0) /
      (G4cbrt(A1[i]) + G4cbrt(AN[i]));
    TM[i] = EEXS - QB - V[i] * A / A1[i];
    if (verboseLevel > 3) {
      G4cout << " i=" << i << " : QB " << QB << " u " << u[i]
         << " V " << V[i] << " TM " << TM[i] << G4endl;
    }
      }
    }
      
    G4double ue = 2.0 * std::sqrt(u1 * EEXS);
    G4double prob_sum = 0.0;

    W[0] = 0.0;
    if (TM[0] > cut_off_energy) {
      G4double AL = getAL(A);
      W[0] = BE * G4cbrt(A1[0]*A1[0]) * G[0] * AL;
      G4double TM1 = 2.0 * std::sqrt(u[0] * TM[0]) - ue;

      if (TM1 > huge_num) TM1 = huge_num;
      else if (TM1 < small) TM1 = small;

      W[0] *= G4Exp(TM1);
      prob_sum += W[0];
    }
      
    for (i = 1; i < 6; i++) {
      W[i] = 0.0;
      if (TM[i] > cut_off_energy) {
    W[i] = BE * G4cbrt(A1[i]*A1[i]) * G[i] * (1.0 + CPA[i]);
    G4double TM1 = 2.0 * std::sqrt(u[i] * TM[i]) - ue;

    if (TM1 > huge_num) TM1 = huge_num;
    else if (TM1 < small) TM1 = small;

    W[i] *= G4Exp(TM1);
    prob_sum += W[i];
      }
    }

    // fisson part
    W[6] = 0.0;
    if (A >= 100.0 && fission_open) {
      G4double X2 = Z * Z / A;
      G4double X1 = 1.0 - 2.0 * Z / A; 
      G4double X = 0.019316 * X2 / (1.0 - 1.79 * X1 * X1);
      G4double EF = EEXS - getQF(X, X2, A, Z, EEXS);
      
      if (EF > 0.0) {
    G4double AF = u1 * getAF(X, A, Z, EEXS);
    G4double TM1 = 2.0 * std::sqrt(AF * EF) - ue;

    if (TM1 > huge_num) TM1 = huge_num;
    else if (TM1 < small) TM1 = small;

    W[6] = BF * G4Exp(TM1);
    if (W[6] > fisssion_cut*W[0]) W[6] = fisssion_cut*W[0];          

    prob_sum += W[6];
      }
    } 

    // again time to decide what next
    if (verboseLevel > 2) {
      G4cout << " Evaporation probabilities: sum = " << prob_sum
         << "\n\t n " << W[0] << " p " << W[1] << " d " << W[2]
         << " He3 " << W[3] << "\n\t t " << W[4] << " alpha " << W[5]
         << " fission " << W[6] << G4endl;
    }

    G4int icase = -1;

    if (prob_sum < prob_cut_off) {      // photon emission chain
      G4double UCR0 = 2.5 + 150.0 / A;
      G4double T00 = 1.0 / (std::sqrt(u1 / UCR0) - 1.25 / UCR0);

      if (verboseLevel > 3)
    G4cout << " UCR0 " << UCR0 << " T00 " << T00 << G4endl;

      G4int itry_gam = 0;
      while (EEXS > cut_off_energy && try_again) {
    itry_gam++;
    G4int itry = 0;
    G4double T04 = 4.0 * T00;
    G4double FMAX;

    if (T04 < EEXS) {
      FMAX = (T04*T04*T04*T04) * G4Exp((EEXS - T04) / T00);
    } else {
      FMAX = EEXS*EEXS*EEXS*EEXS;
    }; 

    if (verboseLevel > 3)
      G4cout << " T04 " << T04 << " FMAX (EEXS^4) " << FMAX << G4endl;

    G4double S(0), X1(0);
    while (itry < itry_max) {
      itry++;
      S = EEXS * inuclRndm();
      X1 = (S*S*S*S) * G4Exp((EEXS - S) / T00);

      if (X1 > FMAX * inuclRndm()) break;
    };

    if (itry == itry_max) {     // Maximum attempts exceeded
      try_again = false;
      break;
    }

    if (verboseLevel > 2) G4cout << " photon escape ?" << G4endl;

    if (S < EEXS) {     // Valid evaporate
      S /= GeV;             // Convert to Bertini units

      G4double pmod = S;
      G4LorentzVector mom = generateWithRandomAngles(pmod, 0.);

      // Push evaporated particle into current rest frame
      mom = toTheNucleiSystemRestFrame.backToTheLab(mom);

      if (verboseLevel > 3) {
        G4cout << " nucleus   px " << PEX.px() << " py " << PEX.py()
           << " pz " << PEX.pz() << " E " << PEX.e() << G4endl
           << " evaporate px " << mom.px() << " py " << mom.py()
           << " pz " << mom.pz() << " E " << mom.e() << G4endl;
      }

      PEX -= mom;           // Remaining four-momentum
      EEXS -= S*GeV;        // New excitation energy (in MeV)

      // NOTE:  In-situ construction will be optimized away (no copying)
      output.addOutgoingParticle(G4InuclElementaryParticle(mom, photon,
                     G4InuclParticle::Equilib));
      
      if (verboseLevel > 3)
        G4cout << output.getOutgoingParticles().back() << G4endl;
      
      ppout += mom;
    } else {
      if (itry_gam == itry_gam_max) try_again = false;
    }
      }     // while (EEXS > cut_off
      try_again = false;
    } else {            // if (prob_sum < prob_cut_off)
      G4double SL = prob_sum * inuclRndm();
      if (verboseLevel > 3) G4cout << " random SL " << SL << G4endl;

      G4double S1 = 0.0;
      for (i = 0; i < 7; i++) { // Select evaporation scenario
    S1 += W[i];
    if (SL <= S1) {
      icase = i;
      break;
    };
      };

      if (icase < 0) continue;  // Failed to choose scenario, try again

      if (icase < 6) { // particle or light nuclei escape
    if (verboseLevel > 2) {
      G4cout << " particle/light-ion escape ("
         << (icase==0 ? "neutron" : icase==1 ? "proton" :
             icase==2 ? "deuteron" : icase==3 ? "He3" :
             icase==4 ? "triton" : icase==5 ? "alpha" : "ERROR!")
         << ")?" << G4endl;
    }

    G4double uc = 2.0 * std::sqrt(u[icase] * TM[icase]);
    G4double ur = (uc > huge_num ? G4Exp(huge_num) : G4Exp(uc));
    G4double d1 = 1.0 / ur;
    G4double d2 = 1.0 / (ur - 1.0);     
    G4int itry1 = 0;
    G4bool bad = true;

    while (itry1 < itry_max && bad) {
      itry1++; 
      G4int itry = 0;
      G4double EPR = -1.0;
      G4double S = 0.0;

      while (itry < itry_max && EPR < 0.0) {
        itry++;
        G4double uu = uc + G4Log((1.0 - d1) * inuclRndm() + d2);
        S = 0.5 * (uc * uc - uu * uu) / u[icase];
        EPR = TM[icase] - S * A / (A - 1.0) + V[icase];
      }; 

      if (verboseLevel > 3) {
        G4cout << " EPR " << EPR << " V " << V[icase]
           << " S " << S << " EEXS " << EEXS << G4endl;
      }

      if (EPR > 0.0 && S > V[icase] && S < EEXS) {  // real escape
        if (verboseLevel > 2)
          G4cout << " escape itry1 " << itry1 << " icase "
             << icase << " S (MeV) " << S << G4endl;

        S /= GeV;               // Convert to Bertini units

        if (icase < 2) {    // particle escape
          G4int ptype = 2 - icase;
          if (verboseLevel > 2)
        G4cout << " particle " << ptype << " escape" << G4endl;
          
          // generate particle momentum
          G4double mass =
        G4InuclElementaryParticle::getParticleMass(ptype);
          
          G4double pmod = std::sqrt((2.0 * mass + S) * S);
          G4LorentzVector mom = generateWithRandomAngles(pmod, mass);
          
          // Push evaporated particle into current rest frame
          mom = toTheNucleiSystemRestFrame.backToTheLab(mom);
          
          if (verboseLevel > 2) {
        G4cout << " nucleus   px " << PEX.px() << " py " << PEX.py()
               << " pz " << PEX.pz() << " E " << PEX.e() << G4endl
               << " evaporate px " << mom.px() << " py " << mom.py()
               << " pz " << mom.pz() << " E " << mom.e() << G4endl;
          }
          
          // New excitation energy depends on residual nuclear state
          G4double mass_new =
        G4InuclNuclei::getNucleiMass(A1[icase],Z1[icase]);
          
          G4double EEXS_new = ((PEX-mom).m() - mass_new)*GeV;
          if (EEXS_new < 0.0) continue; // Sanity check for new nucleus
          
          PEX -= mom;       // Remaining four-momentum
          EEXS = EEXS_new;
          
          A = A1[icase];
          Z = Z1[icase];          
          
          // NOTE:  In-situ construction optimizes away (no copying)
          output.addOutgoingParticle(G4InuclElementaryParticle(mom,
                     ptype, G4InuclParticle::Equilib));
          
          if (verboseLevel > 3)
        G4cout << output.getOutgoingParticles().back() << G4endl;
          
          ppout += mom;
          bad = false;
        } else {    // if (icase < 2)
          if (verboseLevel > 2) {
        G4cout << " nucleus A " << AN[icase] << " Z " << Q[icase]
               << " escape icase " << icase << G4endl;
          }
          
          G4double mass =
        G4InuclNuclei::getNucleiMass(AN[icase],Q[icase]);
          
          // generate particle momentum
          G4double pmod = std::sqrt((2.0 * mass + S) * S);
          G4LorentzVector mom = generateWithRandomAngles(pmod,mass);
          
          // Push evaporated particle into current rest frame
          mom = toTheNucleiSystemRestFrame.backToTheLab(mom);
          
          if (verboseLevel > 2) {
        G4cout << " nucleus   px " << PEX.px() << " py " << PEX.py()
               << " pz " << PEX.pz() << " E " << PEX.e() << G4endl
               << " evaporate px " << mom.px() << " py " << mom.py()
               << " pz " << mom.pz() << " E " << mom.e() << G4endl;
          }
          
          // New excitation energy depends on residual nuclear state
          G4double mass_new =
        G4InuclNuclei::getNucleiMass(A1[icase],Z1[icase]);
          
          G4double EEXS_new = ((PEX-mom).m() - mass_new)*GeV;
          if (EEXS_new < 0.0) continue; // Sanity check for new nucleus
          
          PEX -= mom;       // Remaining four-momentum
          EEXS = EEXS_new;
          
          A = A1[icase];
          Z = Z1[icase];
          
          // NOTE:  In-situ constructor optimizes away (no copying)
          output.addOutgoingNucleus(G4InuclNuclei(mom,
                        AN[icase], Q[icase], 0.*GeV, 
                        G4InuclParticle::Equilib));
          
          if (verboseLevel > 3) 
        G4cout << output.getOutgoingNuclei().back() << G4endl;
          
          ppout += mom;
          bad = false;
        }       // if (icase < 2)
      }     // if (EPR > 0.0 ...
    }       // while (itry1 ...

    if (itry1 == itry_max || bad) try_again = false;
      } else {  // if (icase < 6)
    if (verboseLevel > 2) {
      G4cout << " fission: A " << A << " Z " << Z << " eexs " << EEXS
         << " Wn " << W[0] << " Wf " << W[6] << G4endl;
    }

    // Catch fission output separately for verification
    fission_output.reset();
    theFissioner.deExcite(makeFragment(A,Z,EEXS), fission_output);

    if (fission_output.numberOfFragments() == 2) {      // fission ok
      if (verboseLevel > 2) G4cout << " fission done in eql" << G4endl;

      // Move fission fragments to lab frame for processing
      fission_output.boostToLabFrame(toTheNucleiSystemRestFrame);

      // Now evaporate the fission fragments individually
      this->deExcite(fission_output.getRecoilFragment(0), output);
      this->deExcite(fission_output.getRecoilFragment(1), output);

      validateOutput(target, output);   // Check energy conservation
      return;
    } else { // fission forbidden now
      fission_open = false;
    }
      }     // End of fission case
    }       // if (prob_sum < prob_cut_off)
  }     // while (try_again

  // this time it's final nuclei

  if (itry_global == itry_global_max) {
    if (verboseLevel > 3) {
      G4cout << " ! itry_global " << itry_global_max << G4endl;
    }
  }

  G4LorentzVector pnuc = pin - ppout;

  // NOTE:  In-situ constructor will be optimized away (no copying)
  output.addOutgoingNucleus(G4InuclNuclei(pnuc, A, Z, 0.,
                      G4InuclParticle::Equilib));

  if (verboseLevel > 3) {
    G4cout << " remaining nucleus \n" << output.getOutgoingNuclei().back()
       << G4endl;
  }


  validateOutput(target, output);       // Check energy conservation
  return;
}            

void G4EquilibriumEvaporator::setVerboseLevel ( G4int  verbose )

virtual

Reimplemented from G4CascadeDeexciteBase.

Definition at line 159 of file G4EquilibriumEvaporator.cc.

References G4CascadeDeexciteBase::setVerboseLevel().

Referenced by G4CascadeDeexcitation::setVerboseLevel().

                                                           {
  G4CascadeDeexciteBase::setVerboseLevel(verbose);
  theFissioner.setVerboseLevel(verbose);
  theBigBanger.setVerboseLevel(verbose);
}

---

The documentation for this class was generated from the following files:

• G4EquilibriumEvaporator.hh
• G4EquilibriumEvaporator.cc