#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 G4String &name)

virtual	~G4VCascadeDeexcitation ()

virtual void	collide (G4InuclParticle bullet, G4InuclParticle target, G4CollisionOutput &globalOutput)

Public Member Functions inherited from G4VCascadeCollider
	G4VCascadeCollider (const G4String &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 G4String &name)

Protected Attributes inherited from G4CascadeDeexciteBase
G4CascadeCheckBalance *	balance

G4int	A

G4int	Z

G4LorentzVector	PEX

G4double	EEXS

G4Fragment	aFragment

Protected Attributes inherited from G4VCascadeCollider
G4String	theName

G4int	verboseLevel

Detailed Description

Definition at line 57 of file G4EquilibriumEvaporator.hh.

Constructor & Destructor Documentation

◆ G4EquilibriumEvaporator()

G4EquilibriumEvaporator::G4EquilibriumEvaporator ( )

Definition at line 152 of file G4EquilibriumEvaporator.cc.

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

◆ ~G4EquilibriumEvaporator()

G4EquilibriumEvaporator::~G4EquilibriumEvaporator ( )

virtual

Definition at line 159 of file G4EquilibriumEvaporator.cc.

159{}

Member Function Documentation

◆ deExcite()

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

virtual

Implements G4VCascadeDeexcitation.

Definition at line 170 of file G4EquilibriumEvaporator.cc.

                                                  {
  if (verboseLevel) {
    G4cout << " >>> G4EquilibriumEvaporator::deExcite" << G4endl;
  }
 
  if (verboseLevel>1) G4cout << " evaporating target: \n" << target << G4endl;
 
  // Implementation of the equilibium evaporation according to Dostrovsky
  //   (Phys. Rev. 116 (1959) 683.
  // Note that pairing energy which is of order 1 MeV for odd-even and even-odd
  //   nuclei is not included
 
  // Element in array:  0 : neutron
  //                    1 : proton
  //                    2 : deuteron
  //                    3 : triton
  //                    4 : 3He
  //                    5 : alpha
 
  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 };  // binding energy 
  const G4int AN[6] = { 1, 1, 2, 3, 3, 4 };                     // mass number 
  const G4int Q[6] =  { 0, 1, 1, 1, 2, 2 };                     // charge
  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];
  G4double u[6];     // Level density for each particle type: (Atarg - Aemitted)*parlev  
  G4double V[6];     // Coulomb energy 
  G4double TM[6];    // Maximum possible kinetic energy of emitted particle          
  G4int A1[6];       // A of remnant nucleus
  G4int Z1[6];       // Z of remnant nucleus
 
  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;
    
  /* Loop checking 08.06.2015 MHK */
  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);
      G4double A13 = G4cbrt(A1[0]);
      W[0] = BE * A13*A13 * 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) {
        G4double A13 = G4cbrt(A1[i]);
        W[i] = BE * A13*A13 * 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 = G4double(Z*Z)/G4double(A);
      G4double X1 = 1.0 - 2.0*G4double(Z)/G4double(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*G4UniformRand();
          X1 = (S*S*S*S) * G4Exp((EEXS - S) / T00);
 
          if (X1 > FMAX*G4UniformRand() ) 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*G4UniformRand();
 
      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 Xmax = (std::sqrt(u[icase]*TM[icase] + 0.25) - 0.5)/u[icase];
        G4int itry1 = 0;
        G4bool bad = true;
 
        while (itry1 < itry_max && bad) {
          itry1++; 
          G4int itry = 0;
          G4double S = 0.0;
          G4double X = 0.0;
          G4double Ptest = 0.0;
 
          while (itry < itry_max) {
            itry++;
 
            // Sampling from eq. 17 of Dostrovsky
            X = G4UniformRand()*TM[icase];
            Ptest = (X/Xmax)*G4Exp(-2.*u[icase]*Xmax + 2.*std::sqrt(u[icase]*(TM[icase] - X)));
            if (G4UniformRand() < Ptest) {
              S = X + V[icase];
              break;
            }
          }
 
          if (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;
}            

Referenced by deExcite().

◆ setVerboseLevel()

void G4EquilibriumEvaporator::setVerboseLevel ( G4int verbose )

virtual

Reimplemented from G4CascadeDeexciteBase.

Definition at line 161 of file G4EquilibriumEvaporator.cc.

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

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

Public Member Functions

Additional Inherited Members

Detailed Description

Constructor & Destructor Documentation

◆ G4EquilibriumEvaporator()

◆ ~G4EquilibriumEvaporator()

Member Function Documentation

◆ deExcite()

◆ setVerboseLevel()