G4StatMFChannel Class Reference

#include <G4StatMFChannel.hh>

Data Structures
struct	DeleteFragment

Public Member Functions
G4bool	CheckFragments (void)
void	CreateFragment (G4int A, G4int Z)
	G4StatMFChannel ()
G4FragmentVector *	GetFragments (G4int anA, G4int anZ, G4double T)
G4double	GetFragmentsCoulombEnergy (void)
G4double	GetFragmentsEnergy (G4double T) const
size_t	GetMultiplicity (void)
	~G4StatMFChannel ()

Private Member Functions
void	CoulombImpulse (G4int anA, G4int anZ, G4double T)
void	FragmentsMomenta (G4int NF, G4int idx, G4double T)
	G4StatMFChannel (const G4StatMFChannel &right)
G4bool	operator!= (const G4StatMFChannel &right) const
G4StatMFChannel &	operator= (const G4StatMFChannel &right)
G4bool	operator== (const G4StatMFChannel &right) const
void	PlaceFragments (G4int anA)
G4ThreeVector	RotateMomentum (G4ThreeVector Pa, G4ThreeVector V, G4ThreeVector P)
void	SolveEqOfMotion (G4int anA, G4int anZ, G4double T)

Private Attributes
G4int	_NumOfChargedFragments
G4int	_NumOfNeutralFragments
std::deque< G4StatMFFragment * >	_theFragments

Detailed Description

Definition at line 40 of file G4StatMFChannel.hh.

Constructor & Destructor Documentation

G4StatMFChannel() [1/2]

G4StatMFChannel::G4StatMFChannel

(

)

Definition at line 46 of file G4StatMFChannel.cc.

                                 : 
  _NumOfNeutralFragments(0), 
  _NumOfChargedFragments(0)
{}

~G4StatMFChannel()

G4StatMFChannel::~G4StatMFChannel

(

)

Definition at line 51 of file G4StatMFChannel.cc.

{ 
  if (!_theFragments.empty()) {
    std::for_each(_theFragments.begin(),_theFragments.end(),
          DeleteFragment());
  }
}

References _theFragments.

G4StatMFChannel() [2/2]

G4StatMFChannel::G4StatMFChannel ( const G4StatMFChannel &  right )

private

Member Function Documentation

CheckFragments()

G4bool G4StatMFChannel::CheckFragments

(

void

)

Definition at line 59 of file G4StatMFChannel.cc.

{
  std::deque<G4StatMFFragment*>::iterator i;
  for (i = _theFragments.begin(); 
       i != _theFragments.end(); ++i) 
    {
      G4int A = (*i)->GetA();
      G4int Z = (*i)->GetZ();
      if ( (A > 1 && (Z > A || Z <= 0)) || (A==1 && Z > A) || A <= 0 ) return false;
    }
    return true;
}

References _theFragments, A, and Z.

Referenced by G4StatMF::BreakItUp().

CoulombImpulse()

void G4StatMFChannel::CoulombImpulse ( G4int  anA, G4int  anZ, G4double  T  )

private

Definition at line 136 of file G4StatMFChannel.cc.

{
  // First, we have to place the fragments inside of the original nucleus volume
  PlaceFragments(anA);
    
  // Second, we sample initial charged fragments momenta. There are
  // _NumOfChargedFragments charged fragments and they start at the begining
  // of the vector _theFragments (i.e. 0) 
  FragmentsMomenta(_NumOfChargedFragments, 0, T);
 
  // Third, we have to figure out the asymptotic momenta of charged fragments 
  // For taht we have to solve equations of motion for fragments
  SolveEqOfMotion(anA,anZ,T);
 
  return;
}

References _NumOfChargedFragments, FragmentsMomenta(), PlaceFragments(), and SolveEqOfMotion().

Referenced by GetFragments().

CreateFragment()

void G4StatMFChannel::CreateFragment	(	G4int	A,
		G4int	Z
	)

Definition at line 72 of file G4StatMFChannel.cc.

{
  if (Z <= 0.5) {
    _theFragments.push_back(new G4StatMFFragment(A,Z));
    _NumOfNeutralFragments++;
  } else {
    _theFragments.push_front(new G4StatMFFragment(A,Z));
    _NumOfChargedFragments++;
  }
    
  return;
}

References _NumOfChargedFragments, _NumOfNeutralFragments, _theFragments, A, and Z.

Referenced by G4StatMFMicroCanonical::ChooseAandZ(), G4StatMFMacroCanonical::ChooseZ(), and G4StatMFMicroPartition::ChooseZ().

FragmentsMomenta()

void G4StatMFChannel::FragmentsMomenta ( G4int  NF, G4int  idx, G4double  T  )

private

Definition at line 211 of file G4StatMFChannel.cc.

{
  G4double KinE = 1.5*T*NF; 
  G4ThreeVector p(0.,0.,0.);
    
  if (NF <= 0) return;
  else if (NF == 1) 
    {
      // We have only one fragment to deal with
      p = std::sqrt(2.0*_theFragments[idx]->GetNuclearMass()*KinE)
    *G4RandomDirection();
      _theFragments[idx]->SetMomentum(p);
    } 
  else if (NF == 2) 
    {
      // We have only two fragment to deal with
      G4double M1 = _theFragments[idx]->GetNuclearMass();
      G4double M2 = _theFragments[idx+1]->GetNuclearMass();
      p = std::sqrt(2.0*KinE*(M1*M2)/(M1+M2))*G4RandomDirection();      
      _theFragments[idx]->SetMomentum(p);
      _theFragments[idx+1]->SetMomentum(-p);
    } 
  else 
    {
      // We have more than two fragments
      G4double AvailableE;
      G4int i1,i2;
      G4double SummedE;
      G4ThreeVector SummedP(0.,0.,0.);
      do 
    {
      // Fisrt sample momenta of NF-2 fragments 
      // according to Boltzmann distribution
      AvailableE = 0.0;
      SummedE = 0.0;
      SummedP.setX(0.0);SummedP.setY(0.0);SummedP.setZ(0.0);
      for (G4int i = idx; i < idx+NF-2; ++i) 
        {
          G4double E;
          G4double RandE;
          do 
        {
          E = 9.0*G4UniformRand();
          RandE = std::sqrt(0.5/E)*G4Exp(E-0.5)*G4UniformRand();
        } 
          // Loop checking, 05-Aug-2015, Vladimir Ivanchenko
          while (RandE > 1.0);
          E *= T;
          p = std::sqrt(2.0*E*_theFragments[i]->GetNuclearMass())
        *G4RandomDirection();
          _theFragments[i]->SetMomentum(p);
          SummedE += E;
          SummedP += p;
        }   
      // Calculate momenta of last two fragments in such a way
      // that constraints are satisfied
      i1 = idx+NF-2;  // before last fragment index
      i2 = idx+NF-1;  // last fragment index
      p = -SummedP;
      AvailableE = KinE - SummedE;
      // Available Kinetic Energy should be shared between two last fragments
    } 
      // Loop checking, 05-Aug-2015, Vladimir Ivanchenko
      while (AvailableE <= p.mag2()/(2.0*(_theFragments[i1]->GetNuclearMass()+
                      _theFragments[i2]->GetNuclearMass())));
      G4double H = 1.0 + _theFragments[i2]->GetNuclearMass()
    /_theFragments[i1]->GetNuclearMass();
      G4double CTM12 = H*(1.0 - 2.0*_theFragments[i2]->GetNuclearMass()
              *AvailableE/p.mag2());
      G4double CosTheta1;
      G4double Sign;
 
      if (CTM12 > 1.) {CosTheta1 = 1.;} 
      else {
    do 
      {
        do 
          {
        CosTheta1 = 1.0 - 2.0*G4UniformRand();
          } 
        // Loop checking, 05-Aug-2015, Vladimir Ivanchenko
        while (CosTheta1*CosTheta1 < CTM12);
      }
    // Loop checking, 05-Aug-2015, Vladimir Ivanchenko
    while (CTM12 >= 0.0 && CosTheta1 < 0.0);
      }
 
      if (CTM12 < 0.0) Sign = 1.0;
      else if (G4UniformRand() <= 0.5) Sign = -1.0;
      else Sign = 1.0;      
        
      G4double P1 = (p.mag()*CosTheta1+Sign*std::sqrt(p.mag2()
                              *(CosTheta1*CosTheta1-CTM12)))/H;
      G4double P2 = std::sqrt(P1*P1+p.mag2() - 2.0*P1*p.mag()*CosTheta1);
      G4double Phi = twopi*G4UniformRand();
      G4double SinTheta1 = std::sqrt(1.0 - CosTheta1*CosTheta1);
      G4double CosPhi1 = std::cos(Phi);
      G4double SinPhi1 = std::sin(Phi);
      G4double CosPhi2 = -CosPhi1;
      G4double SinPhi2 = -SinPhi1;
      G4double CosTheta2 = (p.mag2() + P2*P2 - P1*P1)/(2.0*p.mag()*P2);
      G4double SinTheta2 = 0.0;
      if (CosTheta2 > -1.0 && CosTheta2 < 1.0) {
    SinTheta2 = std::sqrt(1.0 - CosTheta2*CosTheta2);
      }
      G4ThreeVector p1(P1*SinTheta1*CosPhi1,P1*SinTheta1*SinPhi1,P1*CosTheta1);
      G4ThreeVector p2(P2*SinTheta2*CosPhi2,P2*SinTheta2*SinPhi2,P2*CosTheta2);
      G4ThreeVector b(1.0,0.0,0.0);
    
      p1 = RotateMomentum(p,b,p1);
      p2 = RotateMomentum(p,b,p2);
    
      SummedP += p1 + p2;
      SummedE += p1.mag2()/(2.0*_theFragments[i1]->GetNuclearMass()) + 
    p2.mag2()/(2.0*_theFragments[i2]->GetNuclearMass());        
        
      _theFragments[i1]->SetMomentum(p1);
      _theFragments[i2]->SetMomentum(p2);
        
    }
  return;
}

References _theFragments, G4Exp(), G4RandomDirection(), G4UniformRand, CLHEP::Hep3Vector::mag(), CLHEP::Hep3Vector::mag2(), P1, P2, RotateMomentum(), CLHEP::Hep3Vector::setX(), CLHEP::Hep3Vector::setY(), CLHEP::Hep3Vector::setZ(), and twopi.

Referenced by CoulombImpulse(), and GetFragments().

GetFragments()

G4FragmentVector * G4StatMFChannel::GetFragments	(	G4int	anA,
		G4int	anZ,
		G4double	T
	)

Definition at line 118 of file G4StatMFChannel.cc.

{
  // calculate momenta of charged fragments  
  CoulombImpulse(anA,anZ,T);
    
  // calculate momenta of neutral fragments
  FragmentsMomenta(_NumOfNeutralFragments, _NumOfChargedFragments, T);
 
  G4FragmentVector * theResult = new G4FragmentVector;
  std::deque<G4StatMFFragment*>::iterator i;
  for (i = _theFragments.begin(); i != _theFragments.end(); ++i)
    theResult->push_back((*i)->GetFragment(T));
 
  return theResult;
}

References _NumOfChargedFragments, _NumOfNeutralFragments, _theFragments, CoulombImpulse(), and FragmentsMomenta().

GetFragmentsCoulombEnergy()

G4double G4StatMFChannel::GetFragmentsCoulombEnergy

(

void

)

Definition at line 88 of file G4StatMFChannel.cc.

{
  G4double Coulomb =
    std::accumulate(_theFragments.begin(),_theFragments.end(),
                    0.0,
                    [](const G4double& running_total,
                       G4StatMFFragment*& fragment)
                    {
                      return running_total + fragment->GetCoulombEnergy();
                    } );
  //      G4double Coulomb = 0.0;
  //      for (unsigned int i = 0;i < _theFragments.size(); i++)
  //    Coulomb += _theFragments[i]->GetCoulombEnergy();
  return Coulomb;
}

References _theFragments.

Referenced by SolveEqOfMotion().

GetFragmentsEnergy()

G4double G4StatMFChannel::GetFragmentsEnergy

(

G4double

T

)

const

Definition at line 104 of file G4StatMFChannel.cc.

{
  G4double Energy = 0.0;
    
  G4double TranslationalEnergy = 1.5*T*_theFragments.size();
 
  std::deque<G4StatMFFragment*>::const_iterator i;
  for (i = _theFragments.begin(); i != _theFragments.end(); ++i)
    {
      Energy += (*i)->GetEnergy(T);
    }
  return Energy + TranslationalEnergy;  
}

References _theFragments.

Referenced by G4StatMF::CalcEnergy().

GetMultiplicity()

size_t G4StatMFChannel::GetMultiplicity ( void  )

inline

Definition at line 64 of file G4StatMFChannel.hh.

{ return _theFragments.size();}

References _theFragments.

Referenced by G4StatMF::BreakItUp().

operator!=()

G4bool G4StatMFChannel::operator!= ( const G4StatMFChannel &  right ) const

private

operator=()

G4StatMFChannel & G4StatMFChannel::operator= ( const G4StatMFChannel &  right )

private

operator==()

G4bool G4StatMFChannel::operator== ( const G4StatMFChannel &  right ) const

private

PlaceFragments()

void G4StatMFChannel::PlaceFragments ( G4int  anA )

private

Definition at line 155 of file G4StatMFChannel.cc.

{
  G4Pow* g4calc = G4Pow::GetInstance();
  const G4double R0 = G4StatMFParameters::Getr0();
  G4double Rsys = 2.0*R0*g4calc->Z13(anA);
 
  G4bool TooMuchIterations;
  do 
    {
      TooMuchIterations = false;
    
      // Sample the position of the first fragment
      G4double R = (Rsys - R0*g4calc->Z13(_theFragments[0]->GetA()))*
    g4calc->A13(G4UniformRand());
      _theFragments[0]->SetPosition(R*G4RandomDirection());
    
    
      // Sample the position of the remaining fragments
      G4bool ThereAreOverlaps = false;
      std::deque<G4StatMFFragment*>::iterator i;
      for (i = _theFragments.begin()+1; i != _theFragments.end(); ++i) 
    {
      G4int counter = 0;
      do 
        {
          R = (Rsys - R0*g4calc->Z13((*i)->GetA()))*g4calc->A13(G4UniformRand());
          (*i)->SetPosition(R*G4RandomDirection());
        
          // Check that there are not overlapping fragments
          std::deque<G4StatMFFragment*>::iterator j;
          for (j = _theFragments.begin(); j != i; ++j) 
        {
          G4ThreeVector FragToFragVector = 
            (*i)->GetPosition() - (*j)->GetPosition();
          G4double Rmin = R0*(g4calc->Z13((*i)->GetA()) +
                      g4calc->Z13((*j)->GetA()));
          if ( (ThereAreOverlaps = (FragToFragVector.mag2() < Rmin*Rmin))) 
            { break; }
        }
          counter++;
          // Loop checking, 05-Aug-2015, Vladimir Ivanchenko
        } while (ThereAreOverlaps && counter < 1000);
        
      if (counter >= 1000) 
        {
          TooMuchIterations = true;
          break;
        } 
    }
      // Loop checking, 05-Aug-2015, Vladimir Ivanchenko
    } while (TooMuchIterations);
    return;
}

References _theFragments, G4Pow::A13(), G4RandomDirection(), G4UniformRand, G4Pow::GetInstance(), G4StatMFParameters::Getr0(), CLHEP::Hep3Vector::mag2(), and G4Pow::Z13().

Referenced by CoulombImpulse().

RotateMomentum()

G4ThreeVector G4StatMFChannel::RotateMomentum ( G4ThreeVector  Pa, G4ThreeVector  V, G4ThreeVector  P  )

private

Definition at line 429 of file G4StatMFChannel.cc.

{
  G4ThreeVector U = Pa.unit();
  
  G4double Alpha1 = U * V;
  
  G4double Alpha2 = std::sqrt(V.mag2() - Alpha1*Alpha1);
 
  G4ThreeVector N = (1./Alpha2)*U.cross(V);
  
  G4ThreeVector RotatedMomentum(
                ( (V.x() - Alpha1*U.x())/Alpha2 ) * P.x() + N.x() * P.y() + U.x() * P.z(),
                ( (V.y() - Alpha1*U.y())/Alpha2 ) * P.x() + N.y() * P.y() + U.y() * P.z(),
                ( (V.z() - Alpha1*U.z())/Alpha2 ) * P.x() + N.z() * P.y() + U.z() * P.z()
                );
  return RotatedMomentum;
}

References CLHEP::Hep3Vector::cross(), CLHEP::Hep3Vector::mag2(), P, CLHEP::Hep3Vector::unit(), CLHEP::Hep3Vector::x(), CLHEP::Hep3Vector::y(), and CLHEP::Hep3Vector::z().

Referenced by FragmentsMomenta().

SolveEqOfMotion()

void G4StatMFChannel::SolveEqOfMotion ( G4int  anA, G4int  anZ, G4double  T  )

private

Definition at line 340 of file G4StatMFChannel.cc.

{
  G4Pow* g4calc = G4Pow::GetInstance();
  G4double CoulombEnergy = 0.6*elm_coupling*anZ*anZ*
    g4calc->A13(1.0+G4StatMFParameters::GetKappaCoulomb())/
    (G4StatMFParameters::Getr0()*g4calc->Z13(anA)) - GetFragmentsCoulombEnergy();
  if (CoulombEnergy <= 0.0) return;
  
  G4int Iterations = 0;
  G4double TimeN = 0.0;
  G4double TimeS = 0.0;
  G4double DeltaTime = 10.0;
  
  G4ThreeVector * Pos = new G4ThreeVector[_NumOfChargedFragments];
  G4ThreeVector * Vel = new G4ThreeVector[_NumOfChargedFragments];
  G4ThreeVector * Accel = new G4ThreeVector[_NumOfChargedFragments];
  
  G4int i;
  for (i = 0; i < _NumOfChargedFragments; i++) 
    {
      Vel[i] = (1.0/(_theFragments[i]->GetNuclearMass()))*
    _theFragments[i]->GetMomentum();
      Pos[i] = _theFragments[i]->GetPosition();
    }
 
  G4ThreeVector distance(0.,0.,0.);
  G4ThreeVector force(0.,0.,0.);
  G4ThreeVector SavedVel(0.,0.,0.);
  do {
    for (i = 0; i < _NumOfChargedFragments; i++) 
      {
    force.set(0.,0.,0.); 
    for (G4int j = 0; j < _NumOfChargedFragments; j++) 
      {
        if (i != j) 
          {
        distance = Pos[i] - Pos[j];
        force += (elm_coupling*_theFragments[i]->GetZ()
              *_theFragments[j]->GetZ()/
              (distance.mag2()*distance.mag()))*distance;
          }
      }
    Accel[i] = (1./(_theFragments[i]->GetNuclearMass()))*force;
      }
 
    TimeN = TimeS + DeltaTime;
    
    for ( i = 0; i < _NumOfChargedFragments; i++) 
      {
    SavedVel = Vel[i];
    Vel[i] += Accel[i]*(TimeN-TimeS);
    Pos[i] += (SavedVel+Vel[i])*(TimeN-TimeS)*0.5;
      }
    TimeS = TimeN;
    
    // Loop checking, 05-Aug-2015, Vladimir Ivanchenko
  } while (Iterations++ < 100);
    
  // Summed fragment kinetic energy
  G4double TotalKineticEnergy = 0.0;
  for (i = 0; i < _NumOfChargedFragments; i++) 
    {
      TotalKineticEnergy += _theFragments[i]->GetNuclearMass()*
    0.5*Vel[i].mag2();
    }
  // Scaling of fragment velocities
  G4double KineticEnergy = 1.5*_theFragments.size()*T;
  G4double Eta = ( CoulombEnergy + KineticEnergy ) / TotalKineticEnergy;
  for (i = 0; i < _NumOfChargedFragments; i++) 
    {
      Vel[i] *= Eta;
    }
  
  // Finally calculate fragments momenta
  for (i = 0; i < _NumOfChargedFragments; i++) 
    {
      _theFragments[i]->SetMomentum(_theFragments[i]->GetNuclearMass()*Vel[i]);
    }
  
  // garbage collection
  delete [] Pos;
  delete [] Vel;
  delete [] Accel;
  
  return;
}

References _NumOfChargedFragments, _theFragments, G4Pow::A13(), source.hepunit::elm_coupling, G4INCL::Eta, GetFragmentsCoulombEnergy(), G4Pow::GetInstance(), G4StatMFParameters::GetKappaCoulomb(), G4StatMFParameters::Getr0(), CLHEP::Hep3Vector::mag(), CLHEP::Hep3Vector::mag2(), CLHEP::Hep3Vector::set(), and G4Pow::Z13().

Referenced by CoulombImpulse().

Field Documentation

_NumOfChargedFragments

G4int G4StatMFChannel::_NumOfChargedFragments

private

Definition at line 99 of file G4StatMFChannel.hh.

Referenced by CoulombImpulse(), CreateFragment(), GetFragments(), and SolveEqOfMotion().

_NumOfNeutralFragments

G4int G4StatMFChannel::_NumOfNeutralFragments

private

Definition at line 97 of file G4StatMFChannel.hh.

Referenced by CreateFragment(), and GetFragments().

_theFragments

std::deque<G4StatMFFragment*> G4StatMFChannel::_theFragments

private

Definition at line 95 of file G4StatMFChannel.hh.

Referenced by CheckFragments(), CreateFragment(), FragmentsMomenta(), GetFragments(), GetFragmentsCoulombEnergy(), GetFragmentsEnergy(), GetMultiplicity(), PlaceFragments(), SolveEqOfMotion(), and ~G4StatMFChannel().

---

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

• source/processes/hadronic/models/de_excitation/multifragmentation/include/G4StatMFChannel.hh
• source/processes/hadronic/models/de_excitation/multifragmentation/src/G4StatMFChannel.cc