#include <G4UnstableFragmentBreakUp.hh>

Inheritance diagram for G4UnstableFragmentBreakUp:

Public Member Functions
	G4UnstableFragmentBreakUp ()

virtual	~G4UnstableFragmentBreakUp ()

virtual G4bool	BreakUpChain (G4FragmentVector , G4Fragment ) final

virtual G4double	GetEmissionProbability (G4Fragment *fragment) final

void	SetVerbose (G4int val)

Public Member Functions inherited from G4VEvaporationChannel
	G4VEvaporationChannel (const G4String &aName="")

virtual	~G4VEvaporationChannel ()

virtual G4double	GetEmissionProbability (G4Fragment *theNucleus)=0

virtual void	Initialise ()

virtual G4double	GetLifeTime (G4Fragment *theNucleus)

virtual G4Fragment *	EmittedFragment (G4Fragment *theNucleus)

virtual G4bool	BreakUpChain (G4FragmentVector theResult, G4Fragment theNucleus)

G4FragmentVector *	BreakUpFragment (G4Fragment *theNucleus)

virtual void	Dump () const

virtual void	SetICM (G4bool)

virtual void	RDMForced (G4bool)

void	SetOPTxs (G4int opt)

void	UseSICB (G4bool use)

Additional Inherited Members
Protected Attributes inherited from G4VEvaporationChannel
G4int	OPTxs

G4bool	useSICB

Detailed Description

Definition at line 56 of file G4UnstableFragmentBreakUp.hh.

Constructor & Destructor Documentation

◆ G4UnstableFragmentBreakUp()

G4UnstableFragmentBreakUp::G4UnstableFragmentBreakUp ( )

explicit

Definition at line 58 of file G4UnstableFragmentBreakUp.cc.

                                                     : fVerbose(1)
{ 
  fLevelData = G4NuclearLevelData::GetInstance();
  for(G4int i=0; i<6; ++i) {
    masses[i] = G4NucleiProperties::GetNuclearMass(Afr[i], Zfr[i]);
  }
}

◆ ~G4UnstableFragmentBreakUp()

G4UnstableFragmentBreakUp::~G4UnstableFragmentBreakUp ( )

virtual

Definition at line 66 of file G4UnstableFragmentBreakUp.cc.

67{}

Member Function Documentation

◆ BreakUpChain()

G4bool G4UnstableFragmentBreakUp::BreakUpChain	(	G4FragmentVector *	results,
		G4Fragment *	nucleus
	)

finalvirtual

Reimplemented from G4VEvaporationChannel.

Definition at line 69 of file G4UnstableFragmentBreakUp.cc.

{
  //G4cout << "G4UnstableFragmentBreakUp::EmittedFragment" << G4endl;
  G4Fragment* frag = nullptr;
  
  G4int Z = nucleus->GetZ_asInt();
  G4int A = nucleus->GetA_asInt();
 
  G4LorentzVector lv = nucleus->GetMomentum();
  G4double time = nucleus->GetCreationTime();
 
  G4double mass1(0.0), mass2(0.0);
 
  // look for the decay channel with normal masses
  // without Coulomb barrier and paring corrections
  // 1 - recoil, 2 - emitted light ion 
  if(fVerbose > 1) {
    G4cout << "#Unstable decay " << " Z= " << Z << " A= " << A 
       << " Eex(MeV)= " << nucleus->GetExcitationEnergy() << G4endl;
  }
  const G4double tolerance = 10*CLHEP::eV;
  const G4double dmlimit   = 0.2*CLHEP::MeV;
  G4double mass = lv.mag();
  G4double exca = -1000.0;
  G4bool isChannel = false;
  G4int idx = -1;
  for(G4int i=0; i<6; ++i) {
    G4int Zres = Z - Zfr[i];
    G4int Ares = A - Afr[i];
    if(Zres >= 0 && Ares >= Zres && Ares >= Afr[i]) {
      if(Ares <= 4) {
    for(G4int j=0; j<6; ++j) {
      if(Zres == Zfr[j] && Ares == Afr[j]) {
        /*
          G4cout << "i= " << i << " j= " << j << " Zres= " << Zres
          << " Ares= " << Ares << " dm= " << mass - masses[i] - masses[j]
          << G4endl;
        */ 
        G4double delm = mass - masses[i] - masses[j];
        if(delm > exca) {
          mass2 = masses[i]; // emitted
          mass1 = masses[j]; // recoil
              exca  = delm; 
          idx = i;
              if(delm > 0.0) { isChannel = true; }
          break;
        }
      }
    }
      }
      if(isChannel) { break; }
      // no simple channel
      G4double mres = G4NucleiProperties::GetNuclearMass(Ares, Zres);
      G4double e = mass - mres - masses[i];
      // select excited state
      const G4LevelManager* lman = fLevelData->GetLevelManager(Zres, Ares);
      if(lman && e >= 0.0) {
    mass2 = masses[i];
    mass1 = mres + e*G4UniformRand();
    idx = i;
    isChannel = true;
    break;
      } 
      // if physical channel is not identified
      // check excitation energy
      if(e > exca) {
    mass2 = masses[i];
    mass1 = mres;
        if(e > 0.0) { mass1 += e; }
    exca  = e;
    idx   = i;
      }
    }
  }
  G4double massmin = mass1 + mass2;
  if(mass < massmin) {
    if(mass + dmlimit < massmin) { return false; }
    if(fVerbose > 1) {
      G4cout << "#Unstable decay correction: Z= " << Z << " A= " << A 
         << " idx= " << idx
         << " deltaM(MeV)= " << mass - massmin
         << G4endl;
    }      
    mass = massmin;
    G4double e = std::max(lv.e(), mass + tolerance);
    G4double mom = std::sqrt((e - mass)*(e + mass));
    G4ThreeVector dir = lv.vect().unit();
    lv.set(dir*mom, e);
  }
 
  // compute energy of light fragment
  G4double e2 = 0.5*((mass - mass1)*(mass + mass1) + mass2*mass2)/mass;
  e2 = std::max(e2, mass2);
  G4double mom = std::sqrt((e2 - mass2)*(e2 + mass2));
 
  // sample decay
  G4ThreeVector bst  = lv.boostVector();
  G4ThreeVector v = G4RandomDirection();
  G4LorentzVector mom2 = G4LorentzVector(v*mom, e2);
  mom2.boost(bst);  
  frag = new G4Fragment(Afr[idx], Zfr[idx], mom2);
  frag->SetCreationTime(time);
  results->push_back(frag);
 
  // residual
  lv -= mom2;
  Z  -= Zfr[idx];
  A  -= Afr[idx];
    
  nucleus->SetZandA_asInt(Z, A);
  nucleus->SetMomentum(lv);
  return true;
}